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Chapter  1 


Introduction 

1.1  Motivation  and  objectives 

The  recent  fusion  of  decades  of  advancements  in  mathematical  models,  numer¬ 
ical  algorithms  and  computer  architecture  marked  the  beginning  of  a  new  era  in  the 
science  of  simulation.  A  simulation  that  would  have  taken  centuries  on  a  1947  Mark 
I  computer  using  Gaussian  elimination  now  takes  only  a  few  seconds  with  parallel 
multigrid  on  the  IBM  Blue  Gene/P  supercomputer.  Gomputational  fluid  dynamics 
(GFD),  which  is  at  the  forefront  of  computational  mechanics,  in  utilizing  large-scale 
computational  resources  to  tackle  problems  of  increasing  size  and  complexity,  has 
also  greatly  benehted  from  the  above  developments.  Nowadays,  even  commercially 
available,  general  purpose,  GFD  solvers  offer  some  form  of  parallel  computing  capa¬ 
bility.  The  gains,  however,  in  our  ability  to  resolve  the  energy  transfer  in  convection 
dominated  flows  are  not  as  impressive.  During  the  past  30  years  direct  numeri¬ 
cal  simulations  (DNS),  which  are  three-dimensional,  time-dependent  computations 
where  all  scales  of  motion  down  to  the  Kolmogorov  scale  are  resolved,  have  only 
seen  moderate  increases  in  the  simulated  Reynolds  numbers.  In  addition,  all  these 
state-of-the-art  computations  are  usually  performed  on  idealized  settings  (i.e.,  ho¬ 
mogeneous  turbulence  [115],  fully  developed  channel  flows  [26]).  Applications  to 
more  realistic  conhgurations  are  limited  (see  for  example  [30]). 
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Many  problems  of  interest  in  physics  and  engineering  include  complex  geome¬ 
tries  with  moving  and/or  deforming  boundaries.  The  latter  is  amongst  the  most 
challenging  problems  in  computational  mechanics:  the  main  difficulty  arises  from 
the  fact  that  the  spatial  domain  occupied  by  the  fluid  changes  with  time  and  the 
location  of  the  boundary  is  usually  an  unknown  itself  that  depends  on  the  fluid  flow 
and  the  motion  and/or  deformation  of  the  body.  There  is  only  a  limited  number  of 
special  cases  where  established  CFD  codes  can  be  directly  applied  to  fluid-structure 
interaction  (FSI)  problems  with  a  relatively  small  overhead.  The  use  of  moving  ref¬ 
erence  frames  [55],  or  coordinate  transformations  [72]  are  characteristic  examples. 
In  more  complex  conhgurations  there  are  two  classes  of  methods  that  can  be  used: 
i)  formulations  utilizing  moving  and/or  deforming  grids  that  continuously  adapt  to 
the  changing  location  of  the  body  (see,  for  example  [98,  46]);  and  ii)  non-boundary 
conforming  methods,  where  the  requirement  for  the  grid  to  conform  to  the  body  is 
relaxed  and  boundary  conditions  are  imposed  by  using  external  forcing  functions,  or 
local  reconstructions  (see  [68]  for  a  recent  review  of  the  different  strategies).  The  lat¬ 
ter  family  of  methods  has  significant  advantages  in  conhgurations  involving  multiple 
bodies  undergoing  large  motions  and/or  deformations,  compared  with  the  former, 
where  the  need  for  constant  deformation/regeneration  of  the  grid  has  an  adverse 
impact  on  the  accuracy  and  efficiency  of  the  huid  solvers.  Applications  of  either  of 
the  above  strategies  in  massively  parallel  environments  are  yet  to  be  reported. 

The  last  few  years  have  seen  a  paradigm  shift  in  high  performance  computing 
hardware  from  machines  with  few  powerful  processors  that  achieved  performance  by 
increasing  the  clock  rate,  to  systems  with  many  relatively  simpler  processors  running 
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at  moderate  clock  speeds.  In  the  new  paradigm,  the  computing  power  of  a  processor 
node  is  increased  by  including  multiple  processing  units  with  shared  memory  on  the 
same  chip,  leading  to  hne-grained  parallelism.  This  presents  a  unique  opportunity 
to  advance  high-hdelity,  eddy-resolving  CFD  solvers,  to  the  next  level  and  develop 
scalable  algorithms  for  challenging  multiphysics,  multiscale  problems  such  as  fluid- 
structure  interactions. 

The  main  objective  of  the  proposed  work  is  the  development  of  scalable  tools 
and  algorithms  applicable  to  fluid-structure  interactions  in  viscous  incompressible 
flows.  To  achieve  this  objective,  the  focus  of  the  dissertation  will  be  on  the  following 
specihc  aims: 

1.  Development  of  a  robust  immersed-boundary  formulation  with  adaptive  mesh 
refinement  (AMR):  immersed  boundary  methods  are  well  suited  for  multibody 
FSI  problems.  Existing  methods  however,  lack  robustness,  and  are  conhned 
to  low  Reynolds  numbers  because  the  global  grid  rehnement,  which  is  usu¬ 
ally  required  to  resolve  the  sharper  velocity  gradients  associated  with  higher 
Reynolds  numbers,  making  the  computations  prohibitively  expensive.  To  over¬ 
come  these  constraints,  the  following  is  to  be  carried  out:  i)  development  of 
a  robust  immersed  boundary  method  based  on  a  moving  least  squares  (MLS) 
formulation.  This  approach  is  applicable  to  multibody  problems  without  spe¬ 
cial  treatments,  ii)  development  of  an  AMR  rehnement  strategy  to  locally 
rehne  the  computational  mesh  in  areas  of  sharp  velocity  gradients. 

2.  Development  of  strategies  to  perform  large-eddy  simulations  (LES)  of  turbu- 
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lent  and  transitional  flows  within  the  AMR  solver:  In  DNS  with  AMR,  the 
solution  is  smooth  on  the  grid  scale  and  the  interpolation  errors  due  to  the 
reconstruction  of  the  fluxes  at  the  non-matching  interfaces  between  coarse  and 
hne  grids  are  small.  In  LES,  on  the  other  hand,  the  flow  held  is  generally  not 
smooth  at  the  smallest  scale  (the  hlter  width  is  not  much  larger  than  the  grid 
size),  and  the  numerical  errors  in  the  interpolation  between  grids  can  be  signif¬ 
icant.  In  addition,  when  a  non-uniform  hlter-width  is  used,  differentiation  and 
hltering  do  not  commute,  and  additional  terms  ( “commutator  errors” )  appear 
in  the  equations  of  motion.  In  the  present  work,  the  author  will  investigate  the 
signihcance  of  the  different  errors  on  the  accuracy  of  the  results  and  develop 
strategies  to  eliminate  them. 

3.  Utilization  of  the  tool  in  applications  to  flapping  flight:  The  developed  tools 
will  be  applied  to  a  set  of  two-dimensional  and  three-dimensional  happing 
hight  problems.  A  very  important  question  nowadays  is  how  does  wing  hex- 
ibility  ahect  the  aerodynamic  characteristics  of  a  wing  for  a  given  prescribed 
kinematics.  The  ehects  of  wing  hexibility  are  evaluated  on  the  performance 
of  a  happing,  hexible  airfoil.  In  three-dimensional  studies  most  of  the  work 
on  happing  hight  today  involves  tethered  insects  with  prescribed  kinematics. 
To  demonstrate  the  capabilities  of  the  proposed  method,  the  author  conducts 
computations  of  insects  in  free  hight  (the  wing  motion  is  prescribed  and  the 
overall  system  ’hies’  in  response  to  the  aerodynamic  forces  that  are  produced). 
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1.2  Prior  work 


The  combination  of  fluid  and  structural  solvers  on  FSI  simulations  requires  a 
strategy  to  impose  boundary  conditions  on  the  fluid,  and  also  dehne  the  fluid  forces 
acting  on  the  solids.  Usually,  the  fluid  is  described  on  an  Eulerian  reference  frame, 
while  it  is  simpler  to  use  a  Lagrangian  description  for  the  immersed  bodies.  A  way  of 
overcomming  this  complication  is  to  use  the  Arbitrary  Lagrangian-Eulerian  (ALE) 
formulation  [29]  for  the  fluid  description.  Here,  the  grid  is  required  to  conform  to 
the  immersed  bodies,  in  the  sense  that,  the  boundaries  of  these  are  also  boundaries 
of  the  discrete  fluid  domain.  Then,  the  vertices  (or  nodes)  of  the  fluid  mesh  are 
displaced  by  using  a  conveniently  dehned  velocity  held,  to  account  for  the  motion  of 
internal  surface  boundaries  of  the  huid  domain.  Boundary-conforming  methods  have 
the  disadvantage  that,  in  large  motion  or  deformation  regimes,  the  grid  distortion 
incurred  has  a  negative  impact  on  the  accuracy  of  calculations  [8].  It  is  possible 
to  combine  ALE  schemes  with  remeshing,  for  cases  where  the  mesh  has  reached  an 
unacceptable  level  of  distortion,  but  at  the  overhead  of  modifying  the  grid  topology 
and  implementing  complicated  interpolation  operations  for  grid  variables. 

On  the  other  hand,  immersed  boundary  methods  can  be  used  to  compute 
the  how  around  immersed  bodies  by  solving  the  Navier-stokes  equations  of  motion 
on  a  hxed  structured  grid,  generally  not  aligned  with  the  body.  Depending  on 
the  specihcs  of  the  formulation,  boundary  conditions  are  imposed  by  appropriately 
modifying  the  stencil  in  the  neighborhood  of  the  body  [100],  or  by  using  a  forc¬ 
ing  function  on  the  Navier-Stokes  equations  which  can  be  derived  either  by  using 
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physical  arguments  [76],  or  directly  from  the  discrete  problem  [37]. 

In  the  case  of  Peskins  [76]  formulation,  the  presence  of  the  body  in  the  flow 
is  taken  into  account  through  an  elastic  deformation  force,  which  is  transferred  to 
the  fluid  momentum  equation  as  a  pseudo  body  force.  This  transfer  requires  special 
treatment  and  it  is  done  by  using  kernel  functions.  A  weakness  of  this  method  is 
that  limiting  inelastic  body  cases  result  in  numerically  stiff  problems,  requiring  very 
small  integration  time  steps.  Several  modihcations  to  the  above  method  have  been 
suggested.  Goldstein  et.  al  [41]  applied  a  forcing  term  governed  by  a  feedback  loop, 
on  the  momentum  equations  discretized  by  a  spectral  method.  They  were  able  to 
effectively  reduce  the  timestep  constraint  inherent  to  Peskins  scheme.  Saiki  and 
Biringen  [86]  extended  the  method  of  Goldstein  et.  al  [41]  to  higher-order  hnite 
difference  schemes  and  successfully  applied  it  to  simulations  of  low  Re  flow  around 
rigid  cylinders.  Peskins  immersed  boundary  method  has  been  applied  to  different 
types  of  low  Reynolds  number  biological  flows  [76],  [77],  [78]. 

A  different  approach  was  taken  by  Glowinski  et  al.  ([39],  [40])  ,  who  imple¬ 
mented  a  fictitious  domain  method,  enforcing  the  fluid  velocity  to  the  given  value 
inside  the  rigid  particle  by  a  distributed  Lagrange  multiplier.  Baaijens  [8]  devel¬ 
oped  a  distributed  Lagrange  multiplier  (DLM),  mortar  element  hctitious  domain 
method  where  the  no  slip  condition  in  the  boundary  of  the  obstacle  is  imposed  as 
an  equation  for  the  Lagrange  multiplier  on  this  boundary.  He  applied  the  method 
to  two-dimensional  fluid-structure  interactions  of  flexible  membranes  subjected  to 
a  fluctuating  channel  flow.  The  DLM  hctitious  domain  method  was  extended  by 
Yu  [116]  for  the  case  of  huid  hexible  body  interactions  where  the  constitutive  equa- 
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tion  used  is  the  one  of  an  incompressible  neo-Hookean  solid.  In  this  formulation  the 
DLM  is  used  to  enforce  the  motion  of  the  fluid  as  the  solid  not  only  on  the  boundary 
but  also  inside  the  deformable  solid  domain. 

Mohd-Yusof  [69]  and  Fadlun  et  al.  [37]  introduced  a  non-boundary  conform¬ 
ing  formulation  based  on  a  direct  forcing  approach.  In  this  case,  the  forcing  is 
done  in  the  discrete  equations  defined  on  the  Eulerian  reference  frame,  and  it  is  de- 
hned  such  that  the  boundary  conditions  are  enforced  on  the  boundary  itself.  These 
methods  -also  called  direct-forcing  methods-  are  particularly  attractive,  especially 
when  combined  with  hnite-difference  or  hnite-volume  formulations,  since  they  can 
be  easily  implemented  in  a  manner  that  does  not  affect  the  efficiency  and  stability  of 
the  solver.  Boundary  motion,  however,  introduces  additional  complications,  and  a 
straightforward  extension  of  the  direct-forcing  formulations  designed  for  stationary 
boundaries  (see  for  example  [37],  [50],  [9])  to  fluid-structure  interaction  problems, 
leads  to  hydrodynamic  forces  that  lack  smoothness  and  are  a  potential  source  of 
instabilities  [101,  113]. 

Yang  and  Balaras  [113]  suggested  that  the  large  fluctuations  of  the  hydrody¬ 
namic  forces  on  moving  immersed  bodies  were  due  to  the  fact  that,  at  any  given 
timestep,  some  of  the  Eulerian  grid  points  in  the  vicinity  of  the  body  will  not  have 
the  correct  velocity,  pressure  or  their  derivatives,  due  to  their  association  with  the 
solid  in  a  previous  timestep.  The  problematic  cells  as  well  as  the  appropriate  treat¬ 
ment  depends  on  the  details  of  the  implementation.  Yang  and  Balaras  [113],  for 
example,  proposed  a  held-extension  procedure,  where  the  solution  is  ’extended’  into 
the  body  in  a  way  that  the  cells  that  emerge  into  the  fluid  have  the  proper  velocity 
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and  pressure  at  later  timesteps.  Mittal  et  al.  [67]  in  their  generalized  ghost-cell  for¬ 
mulation,  assigned  the  proper  values  at  the  problematic  cells  by  interpolating  from 
their  surroundings.  Yang  et.  al  [114]  combined  the  above  method  with  a  predictor- 
corrector  time  integrator  for  structural  equations,  and  applied  it  to  two-dimensional 
rigid  body  motion  problems.  A  similar  scheme  was  applied  by  De  Tullio  and  col¬ 
laborators  [99]  to  DNS  and  FSI  simulation  of  bileaflet  prostatic  heart  valves  under 
pulsatile  flows  and  physiological  conditions. 

Uhlmann  [101]  suggested  an  alternative  direct-forcing  scheme,  where  the  force 
is  computed  on  the  Lagrangian  markers  rather  than  Eulerian  points  as  it  was  done 
in  all  previous  implementations,  which  resulted  in  much  smoother  hydrodynamic 
forces.  He  applied  the  Lagrangian  forcing  scheme  to  low  Re  FSI  particle  sedimenta¬ 
tion  problems.  A  major  drawback  of  immersed  boundary  schemes  coupled  to  a  single 
logically  cartesian  grid,  specially  in  moving  bondaries,  is  that  the  local  resolution 
required  in  areas  where  the  bodies  are  present  is  unavoidably  extended  to  distant 
regions  of  the  domain.  In  order  to  overcome  this  issue  adaptive  mesh  rehnement  of 
the  fluid  grid  can  be  employed. 

Over  the  past  decades  a  signihcant  amount  of  work  on  adaptive  meshing  has 
been  done  in  the  framework  of  unstructured  grid  solution  methods.  The  lack  of 
inherent  structure  of  such  grids  usually  allows  for  a  straightforward  implementation 
of  a  variety  of  AMR  strategies.  A  widely  used  approach  is  the  so  called  h-rehnement, 
where  local  rehnement  is  achieved  by  splitting  existing  cells  into  several  smaller  ones, 
or  by  locally  introducing  additional  nodes  (see  [65]  for  a  review).  For  problems 
that  involve  moving  boundaries,  however,  local  mesh  motion  (r-rehnement)  is  also 


necessary  to  maintain  grid  quality  [98].  A  drawback  of  the  above  methods,  especially 
in  problems  with  large  boundary  motions  and  deformations,  is  the  difficulty  to 
control  grid  quality,  which  has  an  adverse  impact  on  accuracy  and  stability  of  the 
computations. 

Two  general  categories  of  AMR  for  structured  grids,  which  can  be  ideally  cou¬ 
pled  to  the  class  of  non-boundary  conforming  methods  mentioned  above,  can  be 
identihed:  i)  isotropic  splitting  of  individual  cells  that  can  be  managed  using  hi¬ 
erarchical  tree  [15],  or  fully  unstructured  data-structures  [45];  ii)  grid  embedding, 
where  block-structured  grids  composed  of  nested  rectangular  patches  are  used.  The 
latter  approach  maintains  most  of  the  advantages  of  structured  grid  methods  and 
is  an  attractive  platform  for  introducing  AMR  capabilities  in  eddy  resolving  tech¬ 
niques  such  as  the  LES  and  DNS.  Structured  adaptive  mesh  rehnement  (S-AMR) 
was  initially  introduced  by  Berger  and  Oliger  [20]  for  the  solution  of  one-  and  two- 
dimensional  hyperbolic  problems,  where  the  sharp  discontinuities  in  the  solution 
were  better  captured  with  increasingly  rehned  rectangular  grid  patches.  Since  then, 
the  method  has  been  extended  to  three-dimensions,  and  it  has  been  demonstrated 
to  be  a  robust,  cost  effective  approach  for  a  range  of  hyperbolic  problems  (see  for  ex¬ 
ample  [16],  [19]).  Applications  to  incompressible  flows,  however,  have  been  limited, 
primarily  due  to  complications  associated  to  the  enforcement  of  the  divergence-free 
constraint. 

Most  algorithms  for  incompressible  flows  are  based  on  the  extension  of  the 
second-order  projection  method  by  Bell  and  collaborators  [17]  on  the  grid  topology 
proposed  in  [20].  In  this  particular  splitting  scheme  the  viscous  and  advective  terms 
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in  the  momentum  equation  are  advanced  via  a  Crank-Nicolson  scheme,  and  an  un- 
split,  second-order  upwind  Godunov  method  is  used  to  evaluate  the  nonlinear  term 
at  the  time-centered  location.  In  most  cases  time  rehnement  (i.e.  different  time- 
steps  are  used  in  different  AMR  levels)  is  also  employed  and  synchronization  of  the 
solution  between  AMR  levels  is  required  to  ensure  the  divergence  free  constraint. 
Recent  applications  to  multiphase  flows  and  to  prototypical  laminar  flows  can  be 
found  [2]  and  [64]  respectively.  Applications  of  S-AMR  strategies  such  as  the  above 
to  fluid-structure  interaction  problems  have  been  limited  (see  for  example  [42]).  Of 
particular  interest  is  the  work  of  Roma  et.  al  [84] ,  who  applied  the  a  fully  implicit 
version  of  the  immersed  boundary  method  [76]  on  the  two-dimensional  viscous  in¬ 
compressible  flow  equations.  These,  in  turn,  were  solved  using  an  implicit  projection 
method  inspired  by  the  method  given  in  [17]  on  the  composite  grid  structure  of  [20]. 
Their  results  obtained  by  using  self  adaptation  on  the  contraction  of  a  2D  elastic 
spherical  balloon  showed  no  signihcant  difference  between  AMR  calculations  and 
uniform  grid  calculations  with  same  resolution  as  the  hnest  AMR  level. 

Most  of  the  AMR  applications  to  date  for  the  solution  of  turbulent  and  tran¬ 
sitional  flows  have  been  in  DNS  and  Reynolds  Averaged  Navier-Stokes  (RANS)  so¬ 
lutions,  but  their  understanding  and  use  in  LES  of  turbulent  flows  is  reduced.  This 
is  due  to  the  fact  that  in  LES  the  flow  held  is  generally  not  smooth  at  the  smallest 
scale  (the  hlter  width  is  not  much  larger  than  the  grid  size),  so  that  numerical  errors 
in  the  interpolation  between  grids  can  be  signihcant.  In  addition,  the  subgrid-scale 
(SGS)  eddy  viscosity  used  to  represent  the  ehect  of  the  unresolved  scales  is  usu¬ 
ally  proportional  to  the  hlter  width  (and,  in  most  cases,  to  the  grid  size)  squared. 
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Then,  a  sudden  mesh  rehnement  or  coarsening  may  result  in  a  discontinuity  in  eddy 
viscosity,  which  can  generate  significant  errors.  Kravchenko  et  al.  [51],  for  exam¬ 
ple,  computed  plane  channel  flows  using  a  method  based  on  B-splines  that  allowed 
them  to  vary  the  grid  size  (and  hence  the  filter  width)  in  planes  parallel  to  the  wall. 
They  observed  errors  on  the  velocity  field  when  the  coarse-to-fine  grid  ratio  was  2, 
and  noted  that  better  agreement  with  a  single-grid  calculation  could  be  achieved 
either  by  using  a  different  coarse-to-fine  grid-size  ratio,  or  by  providing  a  smoother 
transition  between  meshes.  Of  particular  interest  is  the  work  by  Pantano  et  al.  [75] 
on  a  hybrid  AMR  finite-difference  LES  of  compressible  flows  using  weighted  non- 
oscillatory  schemes.  They  found  that  the  conservation  properties  of  their  fine-coarse 
upwinding  scheme  behave  well  in  cases  where  the  flow  structures  passing  through 
the  interface  are  well  resolved  (he.,  the  flow  is  smooth  on  the  coarse-grid  level). 

Recently,  two  studies  have  focused  on  evaluating  the  errors  due  to  variable 
or  discontinuous  filters,  examining  the  combined  effect  of  commutator  errors  and 
discontinuous  eddy  viscosity.  Cubero  and  Piomelli  [25]  performed  LES  of  plane 
channel  flow  in  which  the  grid  was  uniform  in  the  streamwise  {x)  and  spanwise  {z) 
directions,  and  smoothly  stretched  in  the  wall-normal  direction  y,  however,  a  non- 
uniform  filter  was  imposed.  They  found  that,  as  the  filter  width  is  increased,  the 
wall  stress  and  the  magnitude  of  the  SGS  stress  decreases,  and  the  smaller  scales 
contain  less  energy,  but  the  size  of  the  large  scales  is  essentially  unaffected.  Inverse 
trends  were  observed  for  the  decreasing  filter- width  ratio.  They  also  noted  that,  in 
the  near-wall  region,  the  transition  from  one  filter-width  to  the  other  often  included 
overshoots  or  undershoots.  This  effect  may  be  related  to  the  additional  closure  terms 
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on  the  standard  LES  equations  introduced  by  the  variable  hlter  width;  according  to 
van  der  Bos  and  Geurts  [102],  this  undesirable  phenomenon  may  be  decreased  by 
reducing  the  gradient  of  the  hlter  width. 

Piomelli  et  al.  [80]  applied  the  GDP  hnite-volume  unstructured  code  to  plane 
channel  how  in  which  grid  discontinuities  were  artihcially  introduced,  the  hne-coarse 
grid  interface  being  placed  either  parallel  or  normal  to  the  main  advection  direction. 
Their  hndings  suggest  that,  for  interfaces  parallel  to  the  how  direction,  the  resolved 
stresses  decrease  as  the  grid  is  coarsened,  while  the  subgrid-scale  (SGS)  ones  increase 
proportionately.  If  the  grid  interface  is  placed  close  to  the  buher  layer,  however,  a 
thicker  sublayer  results.  When  the  interface  is  normal  to  the  main  advection,  they 
observed  a  more  complex  behavior:  a  sudden  grid  coarsening  by  a  factor  of  two  in 
each  direction  resulted  in  aliasing  error,  loss  of  phase  information  between  eddies, 
and  decrease  of  the  resolved  Reynolds  stresses.  None  of  the  SGS  models  tested  was 
capable  of  balancing  this  decrease.  A  coarse-to-hne  interface  was  found  to  have  a 
more  benign  character. 

It  is  quite  clear  that  the  grid  discontinuities  that  are  introduced  in  AMR  meth¬ 
ods  may  lead  to  numerical  errors,  in  some  cases  signihcant  ones.  This  is  especially 
true  if  non-dissipative,  centered  schemes  are  used.  In  Ghapter  4,  we  will  study  the 
development  of  turbulence  in  the  vicinity  and  downstream  of  a  rehnement  interface 
evaluating  the  errors  introduced  by  the  grid  discontinuity,  and  the  distance  required 
for  a  return  to  equilibrium.  The  effect  of  different  LES  strategies  on  interpolation 
and  aliasing  errors  around  rehnement  interfaces  will  be  evaluated. 

From  the  above  literature  survey  it  is  clear  that  efficient,  high  hdelity  numerical 
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tools  for  the  solution  of  complex  large-scale  fluid-structure  interaction  problems  on 
massively  parallel  machines  are  yet  to  be  developed.  Such  tools  are  critical  for 
conducting  the  next  generation  engineering  simulations. 

1.3  Outline 

The  remainder  of  this  document  is  organized  as  follows:  In  Chapter  2,  the 
mathematical  model  for  fluid  and  structural  equations  is  laid  out  and  the  developed 
immersed  boundary  scheme  is  explained.  Its  accnracy  and  efficiency  is  evalnated 
using  a  set  of  prototypical  problems  of  increasing  complexity.  In  Chapter  3,  the 
FSI  S-AMR  strategy  is  developed.  The  level  of  accuracy  of  the  required  variable 
interpolation  operators  is  studied,  and  a  novel  divergence-preserving  prolongation 
scheme  for  velocities  is  presented.  Extensive  tests  are  performed  on  a  variety  of 
problems.  In  Chapter  4,  the  effects  of  grid  discontinuities  in  the  case  of  LES  is 
investigated.  A  strategy  based  on  explicit  hltering  of  the  advective  term,  which  is 
effective  in  reducing  numerical  errors  across  the  jump,  is  presented.  To  demonstrate 
the  applicability  of  the  method  in  complex  turbulent  flows  the  case  of  the  flow 
aronnd  a  sphere  at  Re  =  10000  is  presented.  In  Chapter  5,  applications  of  the 
scheme  developed  to  two  and  three  dimensional  problems  relevant  to  flapping  flight 
are  discussed.  Finally,  in  Chapter  6,  conclusions  and  directions  for  future  work  are 
given. 
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Chapter  2 


Immersed  boundary  reconstruction 

In  this  Chapter,  the  author  defines  fiuid  and  structure  models  by  formulating 
the  filtered  Navier-Stokes  and  structural  dynamics  equations.  Then,  the  differences 
between  Eulerian  and  Lagrangian  forcing  as  means  of  imposing  fluid  boundary  con¬ 
ditions  in  immersed  boundary  methods  are  illustrated.  A  thorough  description  of 
the  developed  Moving  Least  Squares  reconstruction  technique  for  Lagrangian  forc¬ 
ing  is  provided,  along  with  an  assessment  of  its  efficiency  and  accuracy  on  selected 
examples. 

2.1  Problem  Formulation 

We  consider  dynamical  systems  consisting  of  fiuid  flow  interacting  with  moving 
and  possibly  deforming  immersed  bodies.  The  flow  is  always  incompressible,  and 
transitional  or  turbulent  flow  patters  may  be  present.  An  example  is  shown  in 
Fig.  2.1,  where  a  set  of  solid  bodies,  $1,  ^2  and  <1>3,  interacts  with  the  flow  within  the 
domain  P  bounded  by  Tq.  In  this  setting,  the  dynamics  of  the  fiuid  and  structures 
are  described  by  different  sets  of  equations,  which  need  to  be  solved  as  a  coupled 
system.  On  the  side  of  the  fluid,  the  LES  modeling  framework  is  adopted,  and  the 
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Figure  2.1:  An  example  of  a  dynamical  system,  where  three  solid  bodies  <Fi,  $2  and 
$3  interact  with  the  flow  in  domain  fl,  with  boundary  F^j. 

spatially-hltered  Navier-Stokes  equations  for  incompressible  flow  are  solved: 
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where  Xj  {i  =  1,  2, 3)  are  the  Cartesian  coordinates,  Uj  are  the  resolved  velocity  com¬ 
ponents  in  the  corresponding  directions,  p  is  the  resolved  pressure,  and  fi  represents 
an  external  body  force  held.  The  equations  are  non-dimensional  and  Re  =  UL/u  is 
the  Reynolds  number  {U  is  a  reference  velocity,  L  a  reference  length  scale  and  u  is 
the  kinematic  viscosity  of  the  huid).  In  LES  the  large  scales  are  resolved  directly  as 
in  a  DNS,  and  all  scales  smaller  than  the  hlter  size,  which  is  usually  proportional 
to  the  local  grid  size,  are  modeled.  In  Eq.  (2.1)  the  effect  of  the  unresolved  scales 
appears  in  the  subgrid  scale  (SGS)  stress  term,  Tij  =  UiUj  —  UiUj,  which  needs  to 
be  parameterized.  We  revisit  LES  in  Chapter  4.  Unless  stated  otherwise,  for  the 
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remainder  of  this  dissertation  we  drop  the  overbars  on  variables,  and  assume  that 
they  are  hltered  variables  in  cases  were  LES  is  used. 

On  the  other  hand,  the  motion  of  a  set  of  rigid  bodies  within  the  fluid  domain, 
O,  is  governed  by  a  set  of  ordinary  differential  equations  (ODEs)  of  the  form  [14]: 
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where  qi  =  [  ^2  •  •  •  gn  ]  ^  ^  vector  containing  the  set  of  n  generalized  coor¬ 

dinates  of  the  structural  system  and  q2  =  qi  is  the  set  of  generalized  velocities.  [I] 
is  the  n  X  n  identity  matrix,  [M(q;^)]  is  the  nonlinear  mass  matrix  and  F(q]^,  q2,  t)  is 
a  vector  containing  damping,  rotation  derived,  elastic  and  externally  applied  forces. 
This  last  set  of  forces  in  our  case,  is  composed  by  the  projection  of  the  gravity  force 
and  fluid  tractions  in  the  direction  of  the  coordinates  g*  {i  =  1, . . . ,  n).  We  note  that 
for  deformable  elastic  bodies,  discretization  of  their  governing  equations  would  also 
lead  to  hrst  order  systems  of  the  form  given  by  Eq.  (2.3).  Therefore,  the  solution 
procedure  presented  here  can  be  extended  to  treat  deformable  bodies. 


2.2  Eulerian  and  Lagrangian  forcing  on  immersed  boundary  methods 

Immersed  boundary  methods  provide  the  means  to  impose  no-slip  boundary 
conditions  on  immersed  surfaces  that  do  not  conform  to  the  Eulerian  fluid  grid.  The 
velocity  held  including  the  boundary  effect  is  reconstructed  around  these  surfaces 
by  applying  a  force  held  on  the  right  hand  side  of  the  discrete  momentum  equations. 
This  force  held,  as  noted  in  Capter  1  can  be  dehned  on  the  Eulerian  huid  grid  points 
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([37],  [9],  [67]),  or  on  the  Lagrangian  marker  points  that  dehne  the  immersed  surfaces 
[101].  For  the  purpose  of  illustrating  the  differences  between  the  two  strategies  let 
us  assume  that  Ui  {i  =  1,2,3)  is  a  discrete  approximation  of  the  velocity  held  and 
write  the  time-discretized  form  of  momentum  equation  as: 


u 


n+l 


Ui 


At 


(2.4) 


where  rhs  contains  all  advective,  diffusive  and,  if  performing  LES,  the  SGS  terms. 
Also,  fi  is  the  direct-forcing  function  which  is  different  from  zero  only  at  the  grid 
points  in  the  vicinity  of  the  immersed  body,  and  n,  n  -|-  1  refer  to  the  current 
and  next  timestep  respectively.  In  the  direct  forcing  scheme  proposed  in  [37]  or 
[9],  for  example,  for  every  point  where  fi  7^  0,  one  can  replace  in  Eq.  (2.4) 
with  the  desired  velocity  uf  (usually  determined  by  means  of  interpolation  from  the 
surrounding  nodes),  and  hnd: 


f, 


n+1/2 


ut 


Ui 


At 


(2.5) 


Substituting  fi  back  into  Eq.  (2.4)  the  proper  boundary  condition,  =  uf  is 
recovered.  In  the  formulation  proposed  by  Uhlmann  [101]  on  the  other  hand,  the 
direct  forcing  function  is  computed  on  each  Lagrangian  marker,  rather  than  on  the 
Eulerian  grid  nodes,  as  follows: 

Tjd  _  Tjn 

^  (2.6) 

*  At  ^  ’ 


The  upper  case  symbols  in  Eq.  (2.6)  denote  the  same  variables  as  in  Eq.  (2.5),  but 

at  the  Lagrangian  points  on  the  immersed  body.  Setting  U  =  17^  +  At 

which  is  practically  the  Lagrangian  counterpart  of  the  predicted  velocity,  Ui,  (see 
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Eq.  2.8),  we  can  rewrite  Eq.  (2.6)  as: 

JJ^n+l|2  _  Ui  ^  ^ 

Uhlmann  [101]  computed  the  volume  force  from  Eq.  (2.7)  and  used  the  regularized 
delta  functions  introduced  in  [76]  as  kernels  in  the  transfer  of  variables  between  the 
Eulerian  and  Lagrangian  grids.  The  overall  implementation  was  tailored  to  model 
suspended  rigid  spherical  particles  in  a  laminar  and  turbulent  flows,  where  it  was 
demonstrated  to  be  very  efficient  and  robust.  Direct  extension  to  more  complex 
fluid-structure  interaction  problems  however,  hinges  upon  the  requirement  to  have 
uniform  elements  on  the  surface  of  the  body,  as  well  as  on  the  fact  that  only  in¬ 
tegral  hydrodynamic  forces  could  be  computed.  In  the  following  sections,  based 
on  the  ideas  presented  in  [101],  we  propose  a  direct-forcing  scheme  that  utilizes  a 
versatile  Moving  Least  Square  (MLS)  approximation  to  build  the  transfer  functions 
between  the  Eulerian  and  Lagrangian  grids,  and  can  be  applied  to  arbitrary  mov¬ 
ing/deforming  bodies.  We  will  also  propose  a  method  to  compute  the  local  traction 
forces.  The  overall  formulation  utilizes  very  compact  stencils  and,  without  com¬ 
promising  accuracy  and  robustness,  gives  results  that  are  identical  to  ‘sharp’  direct 
forcing  methods. 

2.3  Methodologies 

The  proposed  formulation  will  be  discussed  in  the  framework  of  a  hnite- 
difference,  fractional-step,  Navier-Stokes  solver  for  incompressible  flow.  The  ad- 
vective,  diffusive  and  SGS  terms  are  advanced  explicitly  using  an  Adams-Bashforth 
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scheme,  and  all  spatial  derivatives  are  discretized  using  central,  second-order,  finite- 
differences  on  a  staggered  grid.  Details  on  the  basic  solver  together  with  applications 
in  a  variety  of  wall-bounded  and  free-shear  flows  can  be  found  in  [79,  10,  9].  In  the 
following  sections  we  will  focus  on  the  proposed  direct-forcing  scheme  and  compu¬ 
tation  of  the  hydrodynamic  forces. 

2.3.1  MLS  reconstruction 

In  the  framework  of  the  above  mentioned  splitting  scheme  we  first  take  a 
provisional  step  to  compute  the  intermediate  velocities,  Ui,  which  do  not  satisfy  the 
incompressibility  constraint  and  the  boundary  conditions  on  the  immersed  body: 

/\f  c^tP' 

fi.  =  ^  -  At£-,  (2.8) 

where  77  is  a  discrete  operator  representing  the  spatially  discretized  convective, 
viscous  and  SGS  terms.  Next,  we  will  build  a  direct-forcing  function  that  will 
enforce  the  proper  boundary  conditions  on  all  the  Eulerian  grid  nodes  influenced 
by  the  immersed  body.  As  in  [101]  we  will  compute  the  forcing  function  on  the 
Lagrangian  markers  and  then  transfer  it  to  the  Eulerian  grid  nodes.  Our  transfer 
operators,  however,  will  be  constructed  using  MLS  shape  functions  with  compact 
support  [52,  57].  To  facilitate  this  process,  for  each  Lagrangian  marker  we:  i)  Identify 
the  closest  Eulerian  grid  node.  Referring  to  Fig.  2.2a,  for  example,  the  marker  la 
is  associated  to  the  grid  node  which  is  in  the  center  of  a  cell  with  dimensions 
hx  and  hy  in  the  x  and  y  directions  respectively.  Marker  lb  is  associated  to  the  grid 
node  Xb  and  so  on.  Note  that  more  than  one  Lagrangian  markers  from  the  same. 
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or  different  immersed  bodies,  can  be  associated  with  the  same  Eulerian  grid  node, 
ii)  Dehne  a  support-domain  around  each  Lagrangian  marker,  in  which  the  shape 
functions  will  be  constructed.  In  our  case  the  support  domain  is  a  rectangular  box 
of  size  2Hx  x  2Hy  X  2Hz  centered  at  the  location  of  the  marker.  Hy  and 
are  different  for  each  marker  and  are  proportional  to  the  local  Eulerian  grid.  We 
found  Hx  =  l-2hx,  Hy  =  1.2hy  and  =  1.2hz  (see  Fig.  2.2a)  to  be  sufficient  for  all 
cases  considered  in  this  study;  iii)  Associate  a  volume,  AV^  =  {A^  is  the  area 

of  the  body  surface  associated  to  marker  /,  and  is  a  local  thickness  that  depends 
on  local  grid  size  and  will  be  dehned  in  the  following  paragraphs)  to  each  marker 
point.  In  Fig.  2.2a  the  volumes  and  for  markers  and  lb  respectively 

are  shown.  There  is  no  overlapping  between  successive  volumes,  AV\  and  the  sum 
of  all  local  A’'  is  equal  to  the  total  area  of  the  immersed  object  surface. 

We  can  now  dehne  the  transfer  operator  that  will  enable  the  computation  of  Ui 
from  the  corresponding  velocities,  hi,  given  by  Eq.  (2.8).  Using  the  MLS  method, 
Ui  for  each  Lagrangian  marker,  I,  can  be  approximated  in  its  support  domain  as 
follows: 

m 

=  P^(x)a(x),  (2.9) 

i=i 

where  p^(x)  is  the  basis  functions  vector  of  length  m,  a(x)  is  a  vector  of  coeffi¬ 
cients,  and  X  is  the  position  of  the  Lagrangian  marker.  We  found  that  a  linear 
basis,  p^(x)  =  [  1  X  y  lu®  a  cost-efficient  choice  and  would  represent  the  held 
variation  for  all  variables  up  to  the  accuracy  of  our  spatial  discretization  scheme. 
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(a)  (b) 

Figure  2.2:  (a)  Definition  of  the  support-domain  for  two  neighboring  Lagrangian 
markers,  I  a  and  Ib,  which  are  color  coded  for  clarity.  Xa  and  Xb  denote  the  closest 
Eulerian  nodes  to  I  a  and  Ib  respectively.  The  corresponding  volumes  AV  are  also 
shown  (dashed  line),  (b)  The  normal  probe  defined  by  the  Lagrangian  marker  I  and 
point,  e  is  shown  together  with  the  support  domain  used  in  the  MLS  approximation. 

To  obtain  the  coefficient  vector,  a(x),  the  following  weighted  L2-norm  is  defined: 

ne 

J  =  ^  PL  (x  -  x^)  [p^(x^)a(x)  -  u’^]  ^  ,  (2.10) 

k=l 

where  x^  is  the  position  vector  of  the  Eulerian  point  k  in  the  interpolation  stencil, 
u\  is  the  variable  defined  in  Eq.  (2.8)  for  grid  point  fc,  and  PF  (x  —  x^)  is  a  given 
weight  function  that  will  be  defined  below,  ne  is  the  total  number  of  grid  points 
in  the  interpolation  stencil,  which  for  the  linear  basis  function  above,  involves  five 
and  seven  points  in  two-  and  three-dimensions  respectively.  For  simplicity  we  set 
the  closest  point  to  the  Lagrangian  marker  to  be  the  center  point  in  the  stencil. 
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Minimizing  J  with  respect  to  a(x)  leads  to  the  following  set  of  equations: 

A(x)  a(x)  =  B(x)  with, 

ne 

A(x)  =  ^VF(x-x‘)p(x‘)p’-(x'=), 

fc=l 

B(x)  =  [W"(x-xi)p(x^)  ■■■  W"(x-x”")p(x"")],  and 

u-  =  [«•  ■■■  uTY-  (2.11) 

The  size  of  matrix  A(x)  depends  on  the  size  of  the  basis  vector,  p(x),  and  it  is  3  x  3 
in  two-dimensions  and  4  x  4  in  three-dimensions,  while  B(x)  is  of  size  3  x  ne  in 

two-dimensions  or  4  x  ne  in  three-dimensions.  Combining  Eqs.  (2.9)  and  (2.11)  one 

can  write  C*  as  follows: 

ne 

^i(^)  =  4(x)M'  =  ^^(x)u,^  (2.12) 

k=l 

where  $(x)  =  p(x)  A(x)"^  B(x)  is  a  column  vector  with  length  ne,  containing 
the  shape  function  values  for  marker  point  /.  Cubic  splines  are  used  for  the  weight 
functions,  hh(x  —  x^),  above,  which  can  be  written  as: 


2/3  -  47^2  +  47^=" 

for 

Tk  <  0.5 

fT(x  -  x^)  =  < 

4/3  +  4:7%^  -  4/37T^ 

for 

0.5  <  tT  <  1.0 

(2.13) 

0 

for 

tT  >  1.0 

where  1%  =  |x  —  x^|/iJj.  These  functions  are  monotonically  decreasing  and  are 
sufficiently  smooth  in  the  support  domain.  The  resulting  shape  functions  reproduce 
exactly  the  linear  polynomial  contained  in  their  basis  and  possess  the  partition  of 
unity  property  X)r=i '5^4^)  =  1  Also,  the  held  approximation  is  continuous  on 
the  global  domain  as  the  MLS  shape  functions  are  compatible. 
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Eq.  (2.12)  will  give  Ui,  which  can  then  be  substituted  into  Eq.  (2.7)  to  obtain 
the  volume  force  E)  on  all  Lagrangian  markers.  To  transfer  E)  to  the  Eulerian 
points  associated  with  each  marker,  /,  the  same  shape  functions  used  in  interpolation 
procedure  can  be  used  if  properly  scaled  by  a  factor  q,  which  will  be  determined 
later.  In  such  case,  the  final  forces  on  the  Eulerian  grid  nodes  would  be: 

nl 

=  (2.14) 

l=l 

where  ff  is  the  volume  force  in  the  Eulerian  point  k  in  the  direction  i,  0).  is  the  shape 
function  previously  obtained  relating  variables  between  grid  point  k  and  marker  /, 
and  Fj:  is  the  force  in  marker  /.  Also,  nl  is  the  number  of  Lagrangian  markers  which 
are  related  to  the  grid  point  k.  To  properly  rescale  the  shape  functions  we  require 
that  the  total  force  acting  on  the  fluid  is  not  changed  by  the  transfer: 

nte  nil 

=  (2.15) 

k=l  1=1 

where,  AV^  =  {hx  x  hy  x  hz)  is  the  volume  associated  with  the  Eulerian  grid 
point  k,  and  AV^  =  is  the  volume  associated  with  the  marker  /,  with  W  = 
1/3  0L(^a:  +  hy  +  hz).  nte  and  nil  is  the  total  number  of  forced  grid  points,  and 

total  number  of  Lagrangian  markers  respectively.  As  our  surface  is  discretized  using 
triangular  elements,  the  area  for  marker  /  is  obtained  by  a  simple  angle  averaging 
process.  Using  (2.14)  in  (2.15)  and  rearranging  the  sums  in  the  left  hand  side  in 
terms  of  the  total  number  of  markers  we  get 

ntl  ne  ntl  ntl 

E  E  =  E  =  E  hAV  (2.16) 

1=1  k=l  1=1  1=1 
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where  AV^''  is  the  averaged  Eulerian  grid  volume  associated  to  the  Lagrangian 
marker  1.  For  Eq.  (2.16)  to  hold  the  scaling  factor  ci  needs  to  be  set  to: 


Cl  = 


(2.17) 


One  can  also  show  that  the  above  scheme  guarantees  the  equivalence  of  total 
torque  between  the  Eulerian  and  Lagrangian  meshes: 


nte 


ntl 


X  X  FfcAE' 


(2.18) 


k=l 


1=1 


For  simplicity  we  will  provide  a  proof  in  two-dimensions,  but  the  extension  to  three- 
dimensions  is  straightforward.  In  particular,  the  two-dimensional  form  of  Eq.  (2.18) 
can  be  written  as: 


nte 

E  Afi  -  v'A 

k=l 

ntl 

=  E  ix‘F‘  -  y'q)  aV, 

1=1 

(2.19) 

nte 

ntl 

k=l 

=  ^X'F^AE', 

(2.20) 

nte 

E  //t  Ar* 

k=l 

ntl 

=  J2x‘FiAV‘, 

1=1 

(2.21) 

Equation  (2.19)  will  hold  if  both  (2.20)  and  (2.21)  hold.  In  the  following  we  will 
consider  proof  of  Eq.  (2.21),  and  similar  arguments  can  be  used  for  Eq.  (2.20).  For 
each  Lagrangian  marker,  can  be  expressed  in  terms  of  the  shape  functions  as 
follows: 

ne 

=  (2.22) 

k=l 

Substituting  (2.22)  and  (2.14)  into  (2.21)  and  reordering  the  sums  in  the  LHS: 

ntl  ne  ntl  ne 

E  E  =  E  E  4/av,  (2,23) 

l=l  k=l  1=1  k=l 
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Inspection  of  Eq.  (2.23)  confirms  that  the  eqnivalence  of  total  torque  will  be  satisfied, 
if  for  each  Lagrangian  marker,  I,  the  following  is  true: 

ne  ne 

c,  (2.24) 

k=\  k=l 

Given  that  q  =  and  assuming  is  constant  on  the  Eulerian 

stencil,  (2.24)  is  trivially  satished.  In  summary,  the  proposed  transfer  operators, 
conserve  momentum  on  both  uniform  and  stretched  grids.  For  torque  to  be  con¬ 
served  the  cell  volume  across  the  stencil  should  be  kept  constant  for  each  marker. 
This  is  satished  in  case  of  uniform  grids.  In  other  situations,  the  departure  from 
equivalence  for  torque  will  depend  on  the  amount  of  stretching  of  the  grid.  Numeri¬ 
cal  experiments  on  two  dimensional  meshes  showed  that  this  difference  is  small  (less 
than  0.5%)  for  10%  grid  stretching  in  each  direction. 

Using  the  forcing  function  from  Eq.  (2.14),  we  can  now  correct  the  intermediate 
velocity  Ui  to  respect  the  boundary  conditions  on  the  immersed  body:  u*  =  Ui+Atfi. 
The  resulting  approximate  velocity  held,  u*,  which  is  not  divergence  free,  can  be 
projected  into  a  divergence-free  space  by  applying  a  correction  of  the  form: 


“P'  =  <  -  Ai  A  (sp) , 


(2.25) 


where  5p  =  is  the  pressure  correction,  which  satishes  the  following  Poisson 

equation: 


(9^  ((5p)  1  du* 

dxidxi  At  dxi 


(2.26) 


The  velocity  held,  given  by  Eq.  (2.25)  is  divergence-free  and  satishes  the 

boundary  conditions  to  the  order  of  O(At^)  [50]. 
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2.3.2  Calculus  of  surface  forces 

In  non-boundary  conforming  formulations  the  fact  that  the  computational  grid 
and  the  surface  of  the  body  are  almost  never  aligned,  introduces  complications  to 
the  computation  of  hydrodynamic  forces  generated  by  the  surrounding  fluid.  In  the 
present  formulation  for  the  case  of  rigid  bodies  the  distributed  forcing  function  given 
by  Eq.  (2.14)  can  be  utilized  to  compute  the  total  hydrodynamic  force  on  a  solid 
object,  provided  that  all  interior  points  are  properly  treated  (see  for  example  [101]). 
Extension  of  this  approach,  however,  to  the  general  case  of  moving  and/or  deforming 
bodies  is  not  trivial.  In  the  present  study  we  have  developed  a  methodology  where 
the  local  hydrodynamic  force  per  unit  area  on  a  surface  element, 

is  computed  directly  from  the  flow  held  around  the  body.  In  Eq.  (2.27),  is  the 
hydrodynamic  surface  force  in  Xi  direction,  rji  is  stress  tensor,  and  rij  is  the  direction 
cosine  of  the  normal  unit  vector  in  Xj  direction.  The  use  of  Eq.  (2.27)  requires 
knowledge  of  p  and  dui/dxj  on  the  body  surface.  In  the  formulation  outlined  above 
the  boundary  is  dehned  in  a  sharp  manner,  but  the  pressure  and  velocity  helds  are 
forced  to  vary  smoothly  through  the  surface  of  the  body.  Consequently,  the  use  of 
the  same  transfer  functions  to  estimate  p  and  dui/dxj  at  the  Lagrangian  markers 
would  probably  underestimate  the  actual  traction  forces.  This  was  also  verihed  by 
a  series  of  numerical  experiments  we  conducted  for  the  case  of  the  how  around  an 
oscillating  cylinder  below. 

To  avoid  such  problems,  for  each  Lagrangian  marker,  /,  on  the  body  we  create 
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a  normal  probe  by  locating  an  external  point,  e,  at  a  distance,  /i„,  from  the  surface 
(see  Fig.  2.2b).  The  distance  hn  is  proportional  to  the  local  grid  spacing  and  is 
set  to:  hn  =  {hx  +  hy  +  hz)  /3.  To  compute  the  surface  pressure  at  marker,  /,  we 
hrst  compute  the  pressure,  p®,  at  point  e,  using  the  MLS  formulation  described  in 
the  previous  section.  The  support  domain  in  this  case  is  centered  around  point  e 
as  shown  in  Fig.  2.2b.  Next,  the  value  of  dp/dn  is  obtained  from  the  momentum 
equation  normal  to  the  boundary  [113]: 


dp  Du 

dn  Dt 


(2.28) 


where  n  is  the  normal  unit  vector  passing  through  the  marker  /,  and  ^  is  the 
acceleration  of  the  marker.  The  value  of  the  pressure  at  the  surface  is  then  obtained 
from: 

p'  =  p"  -  ^hn  (2.29) 

The  velocity  derivatives,  dUi/dxj,  at  the  location  e  for  each  Lagrangian  marker, 
/,  are  computed  by  differentiating  equation  (2.12): 


dUi  ^  d(j)k 

dxj  ^  dxj^^' 

J  A:=l  J 


(2.30) 


where  d<pk/dxj  comes  from  the  solution  of  an  additional  system  of  equations  similar 
to  (2.11)  [57].  Given  the  fact  that  hn  is  of  the  order  of  the  local  grid  size,  and 
assuming  a  linear  variation  of  the  velocity  near  the  body,  the  derivatives,  dUi/dxj, 
coming  from  equation  (2.30)  are  good  approximation  for  the  derivatives  dui/dxj,  at 
the  surface.  Higher-order  reconstruction  procedures  could  also  be  adopted,  albeit 


at  a  higher  cost.  It  is  shown  next  that  the  above  procedure  reproduces,  the  forces 
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on  the  surface  of  an  immersed  body  very  accurately  when  compared  to  boundary- 
conforming  methods  at  the  same  grid  resolution. 

2.4  Results 

In  this  section,  the  author  presents  two  test  problems  to  demonstrate  the  ac¬ 
curacy  and  robustness  of  the  proposed  formulation.  First,  the  formal  accuracy  is 
examined  for  the  case  of  the  flow  around  a  cylinder  submerged  in  a  driven  cavity. 
Then,  the  flow  around  a  cylinder  oscillating  in  a  cross  flow  is  considered.  Here,  the 
focus  is  on  the  accuracy  of  the  local  force  distribution  on  the  surface  of  the  cylin¬ 
der.  Three  dimensional  examples  including  this  immersed  boundary  formulation  are 
shown  in  the  next  Chapter. 

2.4.1  Accuracy  study 

To  evaluate  the  spatial  accuracy  of  proposed  algorithm  we  performed  simula¬ 
tions  of  the  flow  around  a  cylinder  immersed  in  a  lid-driven  cavity.  In  Figure  2.3(a) 
the  geometry  and  a  typical  vorticity  distribution  are  shown.  For  all  cases  con¬ 
sidered  the  cylinder  diameter  was  set  to  12  =  0.4Lr,  and  the  Reynolds  number. 
Re  =  UiidL^/v  =  1000,  where  Ln  is  the  cavity  size  and  Uud  the  velocity  of  the  top 
boundary.  The  no-slip  conditions  on  the  surface  of  the  cylinder  were  enforced  using 
the  proposed  MLS  reconstruction. 

Although  an  analytical  solution  for  this  problem  is  not  available,  the  overall 
accuracy  of  the  scheme  can  be  evaluated  by  comparing  the  solution  among  grids  at 
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(a) 


(b) 


(c) 


Figure  2.3:  The  flow  around  a  cylinder  immersed  in  a  lid-driven  cavity,  (a)  Compu¬ 
tational  set-up;  (b)  L2  norm  of  the  error  and  (c)  Linf  norm  of  the  error  as  a  function 
of  the  cell  size  Ax.  (A)  u  velocity,  (■)  v  velocity 

different  resolution.  To  facilitate  this  comparison  on  a  staggered  grid  arrangement 
we  considered  meshes  with  36^,  60^,  108^,  180^  and  540^  nodes.  With  this  choice 
the  finest  one  (540^)  is  the  reference  solution,  and  the  average  and  maximum  er¬ 
rors  on  each  of  the  coarser  grids  is  computed  without  the  need  to  interpolate.  In 
Figure  2.3(6),  (c),  the  L2  and  L^nf  norms  of  the  error  are  shown  as  a  function  the 
spatial  resolution.  Both  errors  decrease  with  a  second-order  slope,  indicating  that 
the  second-order  spatial  accuracy  of  the  Cartesian  solver  maintained. 
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2.4.2  Oscillating  cylinder  in  a  cross-flow 


The  ability  of  non-boundary  conforming  methods  to  properly  capture  the  sur¬ 
face  pressure  and  viscous  stress  distribution  is  of  paramount  importance,  especially 
in  fluid-structure  interaction  problems.  To  assess  the  performance  of  the  proposed 
formulation  the  author  considered  the  case  of  a  transversely  oscillating  cylinder  in  a 
cross-flow.  The  dominant  parameters  are  the  Reynolds  number  Re  =  UooD/u  [Uoo 
is  the  inflow  velocity),  the  forcing  frequency,  /e,  and  amplitude,  Oq,  of  the  oscillation. 
When  /e  varies  around  the  natural  shedding  frequency,  /o,  interesting  phenomena 
occur  due  to  the  complex  energy  transfer  between  the  fluid  and  the  body  [43,  44], 
Capturing  the  detailed  flow  physics  for  this  problem  requires  an  accurate  re¬ 
production  of  the  vorticity  dynamics  on  the  surface  of  the  body  and  is  a  stringent 
test  for  non-boundary-conforming  schemes.  The  parametric  space  we  considered  is 
the  one  used  in  the  experiments  by  Gu  et  ah  [43],  the  boundary-conforming  sim¬ 
ulations  of  Guilmineau  and  Queutey  [44],  and  computations  by  Yang  and  Balaras 
[113],  where  an  embedded-boundary  method  with  a  direct  forcing  scheme  is  used. 
The  motion  of  the  cylinder  is  given  by  y(t)  =  aosin(2n fet) .  We  considered  three 
cases  with  Re  =  185,  Oq  =  0.2/1  and  /e//o  =  1-0,  1.1,  1.2  respectively.  For  all  cases 
the  computational  domain  was  set  to  50/1  x  30/1  in  the  streamwise  and  cross-stream 
directions  respectively,  with  the  cylinder  located  at  10/1  from  the  inflow  boundary. 
Free-slip  conditions  are  used  at  the  freestream  boundaries  and  a  convective  condition 
at  the  outflow  boundary  [74],  We  considered  two  grids  with  500  x  450  and  850  x  750 
nodes,  where  the  resulting  cell  size  around  the  cylinder  was  Ax  =  Ay  =  0.008/1 
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and  Ax  =  Ay  =  0.004D  respectively.  A  series  of  tests  for  flow  over  a  stationary 
cylinder  was  first  conducted  to  examine  the  sensitivity  of  the  results  to  the  grid 
resolution.  The  predicted  mean  and  root-mean-square  (rms)  values  of  the  drag  and 
lift  coefficients  on  the  hne  grid  were  Cd  =  1.377,  =  0.0296  and  =  0.461, 

and  the  corresponding  values  on  the  coarser  grid  are  within  1.5%  of  the  above, 
demonstrating  the  grid  independency  of  the  results. 


Figure  2.4:  Drag  and  lift  coefficients  as  a  function  of  time  for  the  case  of  a  cylinder 
oscillating  in  a  cross-flow:  (a)  fe/fo  =  1-0,  and  (b)  /e//o  =  1-2  ( — Cd  and  -  -  Cl). 
(c)  Comparison  of  force  coefficients,  o  Cd,  A  are  the  present  results 

for  the  hne  grid;  — Cd,  -  -  C"^^,  —  from  reference  [44],  and  -|-  x  C£)™*, 
^  (jr^s  [113].  (d)  Phase  angle  between  lift  force  and  vertical  displacement.  A 
are  the  present  results  on  the  hne  grid;  -  -  [44]  and  x  [113]. 
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The  temporal  evolution  of  the  the  lift  and  drag  coefficients  for  the  case  of 
the  oscillating  cylinder  are  shown  in  Figure  2.4a,  b.  It  is  evident  that  the  pro¬ 
posed  scheme  results  in  a  smooth  variation  of  the  force  coefficients  without  special 
treatments.  In  Figure  2.4c  a  comparison  of  Cd,  and  for  the  different 

excitation  frequencies  is  shown  with  the  corresponding  results  in  the  boundary  con¬ 
forming  computations  in  [44]  and  the  computations  by  using  a  direct-forcing  scheme 
in  [113].  A  similar  comparison  for  the  phase  angle  between  the  lift  coefficient  and 
transverse  displacement  of  the  cylinder  is  shown  in  Figure  2Ad.  In  general  the 
agreement  is  excellent.  The  largest  discrepancy  appears  in  C d  and  is  of  the  order 
of  3.5%.  We  should  also  note  that  the  numerical  resolution  around  the  cylinder  in 
our  computations  is  comparable  to  the  one  in  the  reference  computations,  where 
Ax  ~  0.005T). 

An  important  point  in  the  computation  of  the  hydrodynamics  forces,  especially 
in  fluid-structure  interaction  problems,  is  the  consistency  of  the  total  force  and  mo¬ 
ment  exchanged  between  the  fluid  and  solid  systems  (action-reaction).  Uhlmann  [101], 
for  example,  who  utilizes  a  similar  forcing  scheme,  proposes  a  force  computation  ap¬ 
proach  that  results  in  equivalence  of  integral  forces.  For  the  case  of  the  oscillating 
cylinder  we  compared  the  force  coefficients  obtained  by  direct  integration  of  the 
local  stresses  resulting  from  the  normal  probe  approach,  to  ones  obtained  using  the 
approach  in  [101].  The  agreement  is  very  good  and  the  maximum  difference  is  2.3% 
for  the  coarse  grid  (Ax  =  0.008)  and  1.4%  for  the  hne  one  (Ax  =  0.004)  indicating 
the  consistency  of  our  approach. 

While  mean  force  predictions  is  a  good  indicator  of  the  overall  performance 
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(a)  (b) 


Figure  2.5:  Distribution  of  the  pressure  and  skin  friction  coefficients  Cp  and  Cf 
for  the  case  of  a  cylinder  oscillating  in  a  cross-flow.  The  cylinder  is  located  at  the 
extreme  upper  position.  present  results  for  Ax  =  0.008D,  — present  results  for 
Ax  =  0.004D,  o  body-htted  computations  in  [44],  -  -  non  boundary-conforming 
computations  in  [113].  (a)  /e//o  =  1-0,  and  (b)  fe/ fo  =  1-2. 

of  the  method,  they  do  not  necessarily  translate  into  an  accurate  representation 
of  the  local  forces.  In  hgure  2.5  the  distributions  of  pressure  coefficient,  Cp,  and 
the  skin  friction  coefficient,  Cf,  on  the  cylinder’s  surface  are  shown  for  the  time 
instance  corresponding  to  the  extreme  upper  position.  Results  for  both  grids  are 
included  from  our  computations,  and  are  compared  with  the  corresponding  results 
by  Guilmineau  and  Queutey  [44]  and  Yang  and  Balaras  [113].  The  higher  sensitivity 
oi  Cf  to  the  grid  resolution  results  in  slightly  lower  peak  values  on  our  coarse  grid 
computations.  The  results  on  hner  grid  agree  very  well  with  the  reference  data.  Cp 
is  less  sensitive  to  the  grid  resolution  the  results  on  the  different  grids  are  almost 
indistinguishable. 
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Chapter  3 


Adaptive  mesh  rehnement  for  fluid-structure  interaction  problems 

In  this  Chapter,  a  S-AMR  strategy  is  described  for  finding  solutions  of  the 
Navier-Stokes  equations  in  laminar  and  turbulent  incompressible  flows.  Based  on 
the  ideas  introduced  in  [20],  the  single-block  solver  presented  in  the  previous  chapter 
is  employed  on  a  hierarchy  of  sub-grids  with  varying  spatial  resolution.  Each  of  these 
sub-grid  blocks  has  a  structured  Cartesian  topology,  and  it  is  part  of  a  tree  data 
structure  that  covers  the  entire  computational  domain.  One  of  the  main  features 
of  the  present  implementation  is  the  utilization  of  the  Paramesh  toolkit  [62]  to 
keep  track  of  the  grid  hierarchy,  and  perform  the  required  restriction/prolongation 
and  guard-cell  hlling  operations.  To  compute  the  flow  in  complex  geometries,  the 
proposed  ’direct-forcing’  embedded  boundary  method  is  used. 

In  Section  3.1,  a  description  of  the  AMR  grid  topology  is  given  and  the  treat¬ 
ment  of  the  solution  at  block  boundaries  is  discussed.  A  strategy  to  address  issues 
related  to  mass  conservation  at  interfaces  and  grid  adaptation  is  also  proposed. 
In  Section  3.1.5  an  overview  of  the  fluid-structure  interaction  scheme  employed  is 
provided.  In  Section  3.2,  the  following  examples  are  discussed  to  demonstrate  the 
accuracy  and  robustness  of  the  proposed  formulation:  i)  the  Taylor-Green  vortex,  ii) 
three-dimensional  vortex  ring  impinging  on  a  wall  at  Re  570,  hi)  ESI  of  two  dimen¬ 
sional  falling  plates,  and  iv)  ESI  of  a  sphere  bouncing  against  a  wall  at  Re  =  830. 
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(b) 


Figure  3.1:  (a)  Example  grid  hierarchy  and  nested  sub-blocks.  Three  rehnement 
levels,  and  have  been  used,  (b)  Staggered  grid  arrangement  in  cells 

adjacent  to  a  coarse-hne  interface:  pressure  is  located  at  the  center  and  velocities 
at  the  face.  The  velocity  component  in  the  direction  is  not  shown  for  clarity. 

3.1  Adaptive  mesh  refinement 

3.1.1  Grid  topology 

The  computational  grid  consists  of  a  number  of  nested  grid  blocks  with  nx  x 
ny  X  nz  computational  cells  at  different  refinement  levels,  Vl\  /  =  0, 1, ...,  Imax-  The 
coarsest  grid  blocks  at  level  12°  always  cover  the  entire  computational  domain,  and 
local  refinement  is  achieved  by  the  bisection  of  selected  blocks  in  every  coordinate 
direction.  In  this  process  an  arbitrary  block,  b,  at  level  I,  for  example,  will  be  the 
origin  of  eight  children  blocks  at  level  I  +  1  that  occupy  the  same  volume  as  their 
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parent  block.  We  call  leaf  blocks  the  blocks  at  the  highest  level  of  rehnement  present 
on  a  particular  region  of  the  domain  (blocks  that  have  not  been  rehned,  and  therefore 
have  no  children).  Any  leaf  block  at  level  /  shares  a  common  boundary  with  leaf 
blocks  whose  rehnement  level  may  differ  from  I  by  at  most  one  level.  An  example 
is  shown  in  Fig.  3.1a,  where  an  S-shaped  domain  is  discretized  with  an  AMR  grid. 
Level  hlo  covers  the  entire  domain  and  consists  of  6  blocks  (gray).  The  grid  is  locally 
rehned  near  the  center  of  the  domain  by  adding  two  levels  of  rehnement,  (orange) 
and  (yellow).  As  a  result  the  cell  size  in  the  direction  for  a  block  at  the  hnest 
level,  /  =  2,  is  only  a  quarter  of  the  one  at  /  =  0:  =  Aa;°/4,  where  i  =  1,2,  3. 

Note  that  given  the  above  constraints,  Axl  for  computational  cells  on  adjacent  leaf 
blocks  at  diherent  rehnement  levels  will  diher  by  a  factor  of  two.  The  resulting  grid 
structure  is  managed  using  the  octree  data-structure  in  the  Paramesh  toolkit  [62], 
which  enables  a  robust  implementation  and  straightforward  parallelization  of  the 
proposed  algorithm. 

A  staggered  variable  arrangement  is  used  in  each  grid  block  as  shown  in 
Fig.  3.1b,  where  the  velocity  component  in  the  direction  has  been  omitted  for 
clarity.  As  we  did  in  the  previous  Chapter  we  drop  the  overbar  indicating  hltered 
variables,  and  denote  the  pressure  and  velocities  for  block  b  at  level  I  as 
ufjkmi  ■*  =  1,2,3  respectively,  jkm  are  indices  that  identify  a  grid  cell  within  block 
b,  and  i  refers  to  the  orientation  of  the  velocity  component.  Other  variables  such 
as  the  turbulent  viscosity,  ut,  can  be  similarly  dehned.  We  will  also  denote  P{b) 
as  the  parent  block  of  block  b,  and  C{b,i)  as  the  child  block  of  block  b,  where 
i  =  1, ...,  2'^  and  d  is  the  number  of  dimensions  of  the  problem. 
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3.1.2  Prolongation  and  restriction  operators 


Let  us  now  define  the  prolongation  and  restriction  operators  that  are  needed 
to  transfer  variables  between  parent  and  children  blocks.  The  general  form  of  a 
restriction  operation  of  a  variable  0  from  block  C{b,i)  at  level  /  +  1  to  the  parent, 
P{b),  at  level  I  is  given  by: 

E  (3.1) 

p^q^r=il 

where,  ijk  are  the  indexes  of  a  cell  on  P{b)  containing  the  restricted  variable  and 
i'j'k'  are  the  indexes  of  a  cell  in  b.  The  limits  il  and  el  define  the  interpolation 
stencil,  where  the  stencil  size  and  interpolation  coefficients,  depend  on  the 

interpolation  scheme  used.  The  prolongation  of  a  variable  from  block  P{b)  to  b,  on 
the  other  hand,  is  given  by: 

el 

E  (3.2) 

p^q^r=il 

On  a  staggered  grid,  prolongation  and  restriction  operators  are  dehned  separately  for 
the  flow  variables  collocated  at  the  cell  centers  {p’fkmi  faces 

{u^ijkm)-  For  the  former  we  use  a  dimension-by-dimension  interpolation  strategy  that 
utilizes  second-order  Lagrange  polynomials  to  perform  one-dimensional  sweeps  over 
a  three-dimensional  stencil.  In  these  cases  the  weight  factors,  apqr  and  in  Eqs. 
(3.1)  and  (3.2)  are  products  of  the  corresponding  coefficients  in  the  one-dimensional 
Lagrange  polynomials.  For  the  face-centered  components  we  exploit  the  fact  that 
they  are  co-planar,  and  the  above  strategy  is  employed  in  two-dimensions.  Example 
restriction  operations  in  a  two-dimensional  staggered  grid  are  shown  in  Fig.  3.2a. 
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(a)  (b) 


Figure  3.2:  (a)  Interpolation  stencil  utilized  by  two-dimensional  prolongation  op¬ 
erators.  It  contains  nine  points  for  cell-centered  variables  such  as,  p\p  and  three 
points  for  variables  located  at  the  cell  faces,  such  as  uh  and  v\y]  (b)  Variable  ar¬ 
rangement  for  the  construction  of  divergence-preserving  prolongation  operators  in 
two-  dimensions . 

The  pressure,  p\y  located  at  cell-center,  ij,  of  the  parent  block,  is  interpolated  from 
pressures  on  a  3  x  3  stencil  from  the  corresponding  children  blocks.  The  velocities 
u\j  and  vlj,  however,  are  found  using  the  one- dimensional  stencils  (two-dimensional 
in  three  dimensions)  shown  in  the  hgure.  Note  that  in  Fig.  3.2  the  superscript 
indicating  the  block  has  been  dropped  for  simplicity  and  the  parent  and  children 
blocks  are  identihed  by  their  level  superscripts,  I  and  I  +  1  respectively. 

The  operators  given  by  Eqs.  (3.1)  and  (3.2)  satisfy  the  accuracy  requirements 
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but  do  not  guarantee  the  divergence-free  evolution  of  the  velocity  field  through 
the  AMR  grid.  In  the  present  implementation  this  property  is  particularly  desirable 
whenever  prolongations  and  restrictions  are  performed  as  part  of  the  local  grid  adap¬ 
tivity  (local  rehnement  and  de-re£nement)  from  timestep  to  timestep.  Divergence- 
free  restriction  operators  can  be  constructed  in  a  fairly  straightforward  manner.  In 
our  staggered  grid  arrangement,  for  example,  using  face-averaging  for  all  the  velocity 
components  preserves  the  divergence  of  the  velocity  vector.  Divergence-preserving 
prolongation  of  a  vector  held,  on  the  other  hand  is  not  straightforward.  Martin  et 
al.  [64],  in  their  AMR  implementation  for  the  Navier-Stokes  equations  for  incom¬ 
pressible  flow,  use  a  divergence  cleaning  methodology  employing  an  extra  Poisson 
solution.  Balsara  [12]  in  the  framework  of  his  MHD  solver,  proposed  a  divergence- 
free  prolongation  reconstruction  applicable  to  AMR  meshes  with  any  refinement 
ratio,  utilizing  linear  and  quadratic  polynomials.  In  the  present  work  we  have  de¬ 
veloped  divergence-preserving  prolongation  operators  tailored  to  the  specific  AMR 
topology,  where  the  grid  size  between  consecutive  refinement  levels  can  only  differ 
by  a  factor  of  two.  We  will  describe  the  proposed  scheme  in  two- dimensions  for  the 
configuration  shown  in  Fig.  3.2b.  Let  us  define  the  discrete  divergence,  D,  at  an 
arbitrary  cell  ij  of  the  parent  block  as: 


where  A(,  and  Ay  is  the  grid  spacing  on  the  parent  block  in  the  x  and  y  directions 
respectively.  The  above  is  also  the  target  divergence  for  the  reconstructed  velocity 


field  on  the  corresponding  children  blocks.  We  will  first  determine  the  eight  velocity 
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components,  «•+/,+!,  4y-i’  4-\,i'+n  4t/'+i’  ^^ich 

are  located  at  the  faces  of  the  parent  block  as  shown  in  Fig.  3.2b.  This  is  done  using 
one-dimensional,  quadratic,  mass-flux  preserving  interpolations  on  each  face.  For 
example,  on  the  right  face  of  the  coarse  cell  one  can  assume  that  the  velocity  can 
be  approximated  as,  u{y)  =  Uq  -l-  aiy  +  02?/^.  Then,  utilizing  the  the  known  values. 


u\  j  and  on  that  face,  together  with  the  discrete  mass-flux  conservation 


constraint,  u\  j  =  Q ^  the  following  algebraic  system  can  be  assembled: 


jj_i  -  ao  +  ai—  +  a2  — 


1  3A^  'ty 

9  (“it/  +  “!t/+i)  =  «o  +  (A^) 


+  - 1" 


The  solution  of  (3.4)  yields  the  coefficients  ao,  ai,  0:2,  and  the  unknown  velocities 
at  the  midpoints  of  the  face  can  the  be  found  in  a  straightforward  manner.  For 
example,  =  uq  -l-  5aiAy/4  -|-  a2(5Aj,/4)^.  Note  that  the  system  (3.4)  can  be 
written  in  matrix  form,  and  the  inverse  of  the  resulting  3x3  Vandermode  matrix 
can  be  found  analytically.  The  same  approach  is  used  on  all  faces  of  the  parent  cell. 

Next,  the  remaining  four  interior  velocities, 
shown  in  Fig.  3.2b  are  determined.  For  each  of  the  children  blocks  one  can  write  the 
discrete  divergence  equations  and  set  it  equal  to  to  that  from  (3.3).  For  example, 
for  the  block  (F  —  1,/  -|-  1)  in  Fig.  3.2b: 


1/^+1  _  7/^+^ 


77^+1  _  77^+1 


=  Du, 


where  A^+^  and  A^^^  is  the  grid  spacing  on  the  children  block  in  the  x  and  y 
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directions  respectively.  Using  similar  expressions  for  all  four  blocks  we  can  assemble 


a  system  of  four  equations  with  four  unknowns.  It  is  of  the  form,  Au^  =  b,  where 


A  the  4x4  matrix  and  u^,  b  the  unknowns  and  source  vectors  respectively.  A, 
however,  is  a  matrix  of  rank  three  and  cannot  be  inverted.  To  address  this  issue  we 
will  consider  an  additional  independent  equation,  i.e.  express  ji  (or  any  of  the 
internal  velocities)  using  quadratic  interpolation: 


(3.6) 


The  resulting  system  of  equations  can  be  written  as: 
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(3.7) 


where  5  =  (3  =  This  is  now  an  extended  system  of  hve  equa¬ 
tions  and  four  unknowns  which  can  be  written  as:  =  (Ae^Ae)“^Ae^be.  The 

product  (Ae^Ae)“^Ae^  is  the  4x5  left-pseudo-inverse  of  the  system,  which  can  be 
found  analytically.  The  addition  of  Eq.  (3.6)  changes  the  rank  of  the  system  and 
the  above  best  fit  solution  is  actually  the  only  solution  to  the  system.  There  are 
signihcant  differences  between  the  method  above  presented  and  Balsaras  scheme  [12] 
when  a  factor  of  2  is  used  in  his  reconstruction.  Balsaras  method  starts  by  dehning 
a  piecewise-linear  prohle  for  the  divergenceless  vector  held  across  each  face  of  the 
rectangular  region  (the  coarse  cell  in  our  case).  The  slopes  on  this  prohle  are  dehned 
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via  the  minmod  slope-limiter.  In  our  case  we  define  a  parabolic  velocity  profile  that 
maintains  the  mass  fiux  across  the  face  and  is  defined  using  neighboring  coarse  grid 
velocities.  The  Balsara  method  makes  use  of  complete  multidimensional  quadratic 
polynomials  inside  the  region,  with  constraints  such  that,  the  continuous  divergence 
at  each  location  in  the  region  is  zero.  With  these,  the  values  of  the  vector  field  are 
interpolated  at  the  internal  fine-grid  locations.  In  our  case  once  the  coarse  faces 
fine-grid  velocities  are  interpolated,  they  are  used  to  define  one  (or  more  in  three 
dimensions)  internal  fine  grid  velocity,  by  linear  or  quadratic  one  dimensional  inter¬ 
polation.  The  remaining  internal  fine  grid  velocities  are  obtained  using  the  discrete 
divergence  equations  for  each  fine-grid  cell.  In  our  scheme  the  2”^^  order  accuracy 
is  maintained  mainly  be  the  mass-flux  conserving  quadratic  interpolation  done  on 
the  coarse  cell  faces.  The  extension  of  the  above  procedure  to  three-dimensions  is 
straightforward,  and  the  detailed  equations  are  given  in  Appendix  A. 

3.1.3  Treatment  of  the  block  boundaries 

To  facilitate  the  discretization  of  the  equations  of  motion  at  block  boundaries, 
overlapping  between  neighboring  blocks  is  created  by  means  of  two  layers  of  ghost 
cells.  A  two-dimensional  configuration  for  neighboring  blocks  at  refinement  levels  / 
and  /  -|-  1  is  shown  in  Fig.  3.3.  To  assign  values  to  the  ghost  cells  on  the  grid-block 
at  level  /:  i)  the  solution  at  level  /  -|-  1  is  restricted  to  its  parent  grid,  which  has  the 
same  level  of  refinement  as  the  adjacent  coarse  block,  using  Eq.  (3.1);  ii)  filling  of 
the  ghost  cells  is  done  by  simple  ‘injection’  of  the  corresponding  variable  from  the 
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Figure  3.3:  Ghost  cell  configuration  at  the  interface  between  blocks  at  levels  I  and 
Z  + 1.  Filled  symbols  denote  the  centers  of  interior  cells  and  open  symbols  the  centers 
of  the  ghost  cells.  The  grid  nodes  surrounded  by  the  green  dotted  lines  are  the  ones 
on  the  interpolation  stencil  for  the  ghost  cells  identified  by  a  green  circle. 

parent  grid.  To  assign  values  at  ghost  cells  on  the  grid-block  at  level  /  -|-  1  a  similar 
procedure  is  used,  which  utilizes  the  prolongation  operator  given  by  Eq.  (3.2).  In 
Fig.  3.3  the  grid  nodes  in  the  interpolation  stencil  are  shown  for  both  cases. 

3.1.4  Temporal  integration  scheme 

A  standard,  second-order,  fractional-step  method  is  utilized  for  the  temporal 
integration  of  the  governing  equations  [49].  On  each  leaf  block  at  level  Z  =  0, ...,  Imax 
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A 


Interface 


Figure  3.4:  Velocity  components  contributing  to  the  mass  flux  across  a  coarse-flne 
interface. 

the  predicted  velocity,  u\,  can  be  written  as: 


At 

T 


(3.8) 


where  if  is  a  discrete  operator  containing  the  convective,  viscous  and  SGS  terms, 
the  superscript,  u,  refers  to  the  time  level,  and  At  is  the  timestep,  which  is  the  same 
on  all  refinement  levels,  1.  All  spatial  derivatives  are  approximated  using  second- 
order  central  differences  on  a  staggered  grid.  The  predicted  velocity  held,  u\,  which 
is  not  divergence  free,  can  be  projected  into  a  divergence-free  space  by  applying  a 
correction  of  the  form: 


(3.9) 


where  is  the  pressure  correction,  which  satisfies  the  following 


Poisson  equation: 


 1  du\ 


(3.10) 


dxidxi  At  dxi 
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The  solution  of  the  Poisson  equation,  (3.10),  is  done  using  the  multigrid  al¬ 
gorithm  developed  by  Martin  and  Cartwright  [63].  This  method  has  been  designed 
for  block-structured  adaptive  grids  and  maintains  second  order  spatial  accuracy, 
imposing  continuity  of  6p^  and  its  derivatives  across  interface  jumps.  It  uses  a 
residual-correction  formulation  and  employs  red-black  Gauss  Seidel  (RBGS)  point 
relaxation  at  each  level.  In  our  implementation,  in  order  to  increase  the  performance 
of  the  solver,  we  exploit  the  uniformity  of  the  mesh  at  the  coarsest  level  of  the  V- 
cycle,  Isoive,  which  covers  the  entire  computational  domain.  In  particular,  the  author 
utilizes  a  direct  solver,  rather  than  RBGS  iterations,  to  solve  the  residual-correction 
equation  at  that  level.  The  only  drawback  of  this  strategy  is  that  the  computational 
domain  has  to  be  rectangular  since  the  direct  solver  utilizes  FFT  or  Gosine  trans¬ 
forms,  depending  on  the  boundary  conditions.  This  strategy  was  found  to  reduce 
the  number  of  levels  on  a  V-Gycle  and  relaxation-communication  operations,  while 
maintaining  the  convergence  rate  of  the  multigrid  scheme. 

An  important  issue  with  the  application  of  projection  schemes  to  S-AMR  grids 
is  related  to  conservation  of  mass  at  coarse-hne  block  interfaces.  In  Fig.  3.4  a  two- 
dimensional  example  of  the  interface  between  cells  at  levels  /  and  /  -|-  1  is  shown. 
The  velocity  component,  nh,  normal  to  the  interface  at  the  coarse  level  /,  as  well 
as  the  corresponding  velocity  components,  and  at  hne  level  I  +  1,  are 

unknowns  to  be  determined  during  the  solution  process.  They  need,  however,  to 
satisfy  an  additional  constraint  coming  form  the  conservation  of  the  mass  flux  across 
the  interface: 


u. 


=  n'+My'+i  + 


rj' 


U 


l-\-l 


Ay' 


i  +  1 


(3.11) 
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where  A|/^  is  the  grid  spacing  at  level  I  and  Ay''~^^=Ay’' /2,  is  the  grid  spacing  at 
level  /  +  1.  In  general,  Eq.  (3.11)  will  not  hold  after  the  computation  of  the  veloc¬ 
ities  normal  to  the  rehnement  jump,  resulting  to  a  flux  mismatch  localized  in  grid 
rehnement  interfaces.  To  alleviate  this  problem  one  can  assume  that  the  velocity 
components  on  the  hue  grid  at  level,  I  + 1,  are  more  accurate  than  the  corresponding 
coarse  grid  component,  u\j,  which  can  then  be  corrected  to  satisfy  Eq.  (3.11)  (see  for 
example  [19],  [2]).  In  the  present  formulation  the  flux  correction  is  realized  in  the 
following  manner:  i)  the  intermediate  velocities  at  the  hne  level,  u^^-}  are  computed 
from  Eq.  (3.8);  ii)  the  intermediate  velocity  on  the  coarse  level,  u\j,  is  then  corrected 
using  Eq.  (3.11);  and  hi)  during  the  iterative  solution  of  the  Poisson  equation  (3.10) 
the  gradient  of  6p  normal  to  the  refinement  interface  at  the  u\j  location  is  forced  to 
satisfy  the  following  equation: 


Ay^  = 


(3.12) 


iv)  it  is  now  trivial  to  show  that  the  corrected  velocity  from  Eq.  (3.9)  will  also 


conserve  the  mass  flux  across  the  interface. 


3.1.5  Fluid-structure  interaction  algorithm 

A  fundamental  complication  with  two-way,  fluid-structure  interaction  (ESI) 
problems,  is  that  the  prediction  of  the  flow  and  the  corresponding  hydrodynamic 
loads  requires  knowledge  of  the  motion  of  the  structure  and  vice-versa.  In  the  present 
study  a  strong-coupling  scheme  is  adopted,  where  the  fluid  and  the  structure  are 
treated  as  elements  of  a  single  dynamical  system,  and  all  of  the  governing  equations 
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are  integrated  simultaneously  and  iteratively  in  the  time  domain.  The  method  is 
based  on  Hamming’s  4*^-order,  predictor-corrector  formulation,  which  avoids  the 
evaluation  of  the  hydrodynamic  loads  at  fractional  time  steps.  The  details  of  the 
overall  approach  and  a  demonstration  of  the  accuracy  and  efficiency  for  a  variety 
of  fluid-structure  interaction  problems  in  viscous  incompressible  flows  can  be  found 
in  [114].  In  Fig.  3.5,  a  flowchart  of  the  overall  algorithm  as  adapted  to  our  AMR 
implementation  is  shown: 

i)  At  the  beginning  of  each  timestep,  and  only  for  the  case  of  LES,  the  SGS 
eddy  viscosity  by  means  of  a  LES  model.  LES  will  be  discussed  in  next  Chapter. 

ii)  Compute  the  provisional  velocity,  u\,  which  does  not  satisfy  boundary  con¬ 
ditions  on  the  immersed  boundary,  and  is  not  divergence  free.  Assign  values  for  u\ 
to  the  ghost  cells  as  described  in  Section  3.1.3. 

iii)  Perform  AMR  if  necessary  (every  n  steps).  In  particular,  all  leaf-blocks  are 

examined  and  flagged  for  rehnement  or  derehnement,  according  to  specihed  criteria. 
Eor  example,  in  all  ESI  cases  reported  below,  a  leaf-block  is  selected  for  rehnement 
when  it  contains  an  immersed  body,  or  when  the  vorticity  magnitude  is  larger  than 
a  predetermined  threshold.  A  leaf-block  is  marked  for  derehnement,  on  the  other 
hand,  when  an  immersed  body  is  not  present  within  the  block,  and  the  vorticity  is 
below  a  threshold.  Other  criteria,  such  as  velocity  error  norms,  SCS  dissipation  etc., 
could  also  be  used  for  this  purpose.  In  the  case  of  rehnement,  the  newly  created 
children  blocks  are  added  to  the  octree.  All  variables  other  than  u\  and  the  discrete 
operator  in  (3.8),  are  interpolated  using  Eq.  (3.2).  For  u\  and 

the  divergence  preserving  prolongation  procedure  discussed  in  Section  3.1.2  is  used. 
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Figure  3.5:  Summary  of  fluid-structure  interaction  strategy. 

In  this  manner,  the  intermediate  velocities  for  the  time-step  following  the 
remeshing  step  maintain  the  required  divergence  level.  For  each  leaf-block  flagged  for 
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derefinement,  we  first  check  if  its  parent  block  contains  children  that  are  all  flagged 
for  derefinement.  In  such  case  a  restriction  step  is  performed  and  the  corresponding 
leaf-blocks  are  removed  from  the  octree.  Quadratic  restriction  using  Eq.  (3.1)  is 
utilized  for  all  variables  except  u\  and  where  linear  restriction  is  used. 

The  latter  is  a  divergence  preserving  operation. 

iv)  In  the  case  of  ESI  a  set  of  predictor-corrector  sub-iterations  is  appied, 
as  shown  in  the  figure.  Convergence  is  assumed  when  the  L2  error  on  structures 
generalized  coordinates  and  velocities  is  less  than  a  certain  tolerance.  In  all  ESI 
computations  of  the  present  study,  the  tolerance  is  set  to  10“®. 

3.2  Numerical  studies 

In  this  section,  a  series  of  test  problems  with  increasing  complexity  are  studied 
to  evaluate  the  accuracy  and  robustness  of  the  proposed  AMR  formulation.  Initially 
the  spatial  and  temporal  accuracy  of  the  method  is  demonstrated  for  the  Taylor- 
Green  vortex  problem.  Then,  the  case  of  a  three-dimensional  vortex  ring  impinging 
on  a  wall  is  considered.  Finally  two  cases,  where  the  accuracy  and  efficiency  of  the 
method  in  the  presence  of  immersed  bodies  is  examined,  are  presented;  the  ESI 
problem  of  two  falling  plates,  and  sphere  bouncing  against  a  wall  at  Re  =  830. 

3.2.1  Taylor-Green  vortex 

To  investigate  the  numerical  accuracy  of  the  method  the  Taylor-Green  vortex 
problem  is  considered.  The  flow-held  is  represented  by  an  array  of  periodic  counter- 
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Figure  3.6:  AMR  grid  configuration  with  two  refinement  levels  for  the  case  of  the 
Taylor-Green  vortex.  Pressure  isolines  ranging  from  p  =  —0.4  to  0.4  at  t  =  0.03  are 
also  shown. 

rotating  vortices  that  decay  in  time,  and  has  an  analytical  solution  of  the  form: 

Ua  = cos xsiny,  (3.13) 

Va  =  sin X  cosy,  (3-14) 

Pa  = ^ — (cos  2a; -|- COS  2?/) ,  (3.15) 

where  Ua  and  Va  are  the  velocity  components  in  the  x  and  y  directions  respectively, 
Pa  is  the  pressure,  v  is  the  kinematic  viscosity,  and  t  is  the  time.  The  size  of  our 
computational  domain  was  27r  x  2tt,  requiring  the  use  of  a  mix  of  homogeneous 
Dirichlet  and  Neumann  boundary  conditions  for  the  velocity  components.  Grids 
with  two  levels  of  rehnement,  as  well  as  uniform  ones  are  considered.  In  the  latter 
case  32^,  64^,  128^  and  256^  grid  points  are  used  and  the  ghost-cell  Filing  is  per- 
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formed  by  simple  injection  of  the  data  between  blocks  sharing  a  common  boundary. 
In  the  former  case  the  base  grids  utilize  16^,  32^,  64^  and  128^  computational  points 
on  their  coarse  levels,  and  the  rehnement  is  performed  in  the  second  and  fourth 
quadrants  of  the  domain  (see  Fig.  3.6).  To  investigate  the  impact  of  the  ghost-cell 
hlling  scheme  on  the  overall  accuracy  of  the  method  both  linear  and  quadratic  inter¬ 
polation  schemes  are  considered.  In  all  cases  the  equations  of  motion  are  integrated 
for  a  total  of  T  =  0.03  time  units  using  a  timestep  of  At  =  5.0  x  10“^.  The  L2  norm 
of  the  residual  for  the  Poisson  solution  was  kept  in  the  order  of  10“^^. 

In  Fig.  3.6,  pressure  isolines  at  t  =  0.03  are  shown  for  the  AMR  grid  with 
two  levels  of  rehnement  (level  0  being  32^).  It  is  evident  that  the  main  features  of 
the  how  are  captured.  In  Figs.  3.7(a)  and  (c)  the  Ljn/  error  norm  for  both  velocity 
components  at  t  =  0.03,  is  shown  for  all  cases  as  a  function  of  the  grid  size,  A.  Note 
that  in  all  the  two-level  calculations.  A,  is  computed  from  the  highest  rehnement 
level.  As  expected,  in  all  uniform  grid  cases  second-order  accuracy  is  observed.  In 
the  case  of  AMR  the  use  of  linear  interpolation  for  the  ghost-cell  hlling  reduces  the 
spatial  accuracy,  which  is  now  closer  to  hrst  order.  The  higher  order  interpolation 
scheme  described  in  section  3.1.3,  on  the  other  hand,  maintains  the  second  order  for 
both  velocity  components.  It  is  seen  that  the  maximum  error  incurred  by  applying 
quadratic  interpolation  in  the  hue  grid  region  of  the  two-level  calculations  is  of 
the  order  of  what  is  obtained  for  the  coarse  grid  uniform  calculation.  This  is  due 
to  contamination  of  hne  grid  with  coarse-grid  error  levels  and  also  the  inaccuracy 
derived  from  the  quadratic  guard-cell  hlling  itself. 


51 


(c)  (d) 

Figure  3.7:  Accuracy  study  on  Taylor-Green  vortex  problem:  1®*  and  2”*^  order  slopes 
are  added  for  clarity,  (a)  &  (c)  ||u  — Ualloo  and  ||n  — Ualloo  respectively  as  a  function  of 
grid  spacing  at  f  =  0.03.  -|-  uniform  grid;  A  AMR  grid,  linear  interpolation;  □  AMR 
grid,  quadratic  interpolation,  (b)  &  (d)  ||n  —  Ua\\oo  and  ||p  —  Pa\\oo  as  a  function 
of  the  grid  spacing  at  t  =  0.3  for  the  case  with  dynamic  refinement /derehnement 
performed  every  10  timesteps;  -|-  linear  prolongation;  A  quadratic  prolongation; 
□  divergence-preserving  prolongation. 
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In  all  the  above  computations,  the  grid  is  locally  rehned,  but  does  not  evolve 
with  time,  as  required  in  most  moving  boundary  problems.  To  examine  the  effects 
of  the  dynamic  adaptation  of  the  grid  on  the  accuracy  of  the  solver,  the  author  also 
utilizes  the  Taylor-Green  vortex  problem  with  the  same  setup  as  above.  In  this  case, 
however,  starting  from  a  single  level  grid  we  rehne  every  10  iterations  by  adding  one 
rehnement  level  (Fig.  3.6),  and  then  derehne  every  10  iterations.  Refinement,  for 
example,  is  performed  at  iteration  count  n  =  10,  30,  50, ..,  90,  and  derehnement  at 
n  =  20, 40,  ..,80.  At  iteration  number  100  the  solution  is  compared  to  the  analytical 
one,  and  error  norms  for  the  velocity  and  pressure  helds  are  computed. 

An  important  part  of  the  dynamic  grid  adaptation  is  the  prolongation  of  the 
solution  at  the  newly  generated  children  blocks  in  refinement  areas,  and/or  the  re¬ 
striction  of  the  solution  from  eliminated  leaf-blocks  to  its  parent,  in  derehnement 
areas.  As  we  discussed  in  section  3.1.5,  AMR  is  performed  after  the  computation  of 
the  provisional  velocity,  u\.  One  can  show  that  u\  is  within  O(At^)  of  a  divergence- 
free  held  (see  for  example  [49])  and,  therefore,  its  prolongation/restriction  after  a 
grid  adaptation  step  should  maintain  this  property,  so  the  overall  temporal  accu¬ 
racy  is  not  ahected.  As  we  already  discussed  in  Section  3.1.2,  not  all  prolonga¬ 
tion/restriction  operators  are  divergence  preserving. 

To  illuminate  their  ehects  on  the  accuracy  of  the  solution  we  conducted  compu¬ 
tations  with  three  diherent  prolongation  operators  for  the  velocity  held;  linear  and 
quadratic,  which  do  not  preserve  the  divergence  of  u\,  and  the  quadratic  divergence¬ 
preserving  operator  discussed  in  section  3.1.2.  In  Figs.  3.7(b)  and  (d)  the  Linf  error 
norms  for  the  velocities  and  pressure  are  shown  as  a  function  of  the  grid  size  for  all 
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cases.  It  is  evident  that  the  spatial  accuracy  is  not  affected  by  the  choice  of  prolon¬ 
gation  operator,  and  in  all  cases  second  order  accuracy  maintained.  As  expected, 
the  error  using  linear  prolongation  is  about  an  order  of  magnitude  higher  than  using 
the  quadratic  interpolation  variants. 


(a)  (b) 

Figure  3.8:  Accuracy  study  on  the  Taylor-Green  vortex  problem:  (a)  &  (b)  Variation 
of  ||m  —  Malloo,  and  ||p  —  Palloo  respectively  as  a  function  of  time  for  the  case  with 

dynamic  refinement /derefinement  performed  every  10  timesteps;.  - {black)  linear 

prolongation; - {red)  quadratic  prolongation;  .  {blue)  divergence-preserving 

prolongation. 


The  evolution  of  the  error  norms  with  respect  to  time,  on  the  other  hand,  re¬ 
veals  some  very  interesting  patterns,  as  seen  in  Figs.  3.8(a)  and  (b).  The  prolonga¬ 
tion  operators,  which  do  not  preserve  the  divergence  of  u[,  result  in  a  monotonically 
increasing  error  for  both  velocity  components,  while  for  the  divergence-preserving 
prolongation  operator  the  error  is  always  0(10“^).  These  effects  are  more  profound 
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in  the  behavior  of  the  pressure,  where  the  linear  and  quadratic  prolongations  result 
in  an  error  jump  of  at  least  three  orders  of  magnitude  just  after  each  rehnement 
step,  compared  to  the  divergence  preserving  scheme.  The  pressure  recovers  in  a 
few  timesteps,  but  overall  the  error  grows  over  time.  In  the  case  of  fluid-structure 
interaction  problems  these  pressure  oscillations  can  introduce  spurious  loading  on 
the  structure  leading  to  large  errors  on  the  solution.  For  the  divergence-preserving 
scheme  the  error  jump  is  very  small,  indicating  that  the  spurious  oscillations  are 
practically  eliminated.  It  is  important  to  note  that,  the  choice  of  other  projection 
scheme  variants,  where  the  pressure  is  not  treated  incrementally,  or  making  use  of 
high-resolution  approximate  projections  [3]  might  have  a  beneficial  effect  in  damp¬ 
ing  the  observed  pressure  oscillations  and  numerical  errors  due  to  non-divergence 
preserving  interpolations. 

3.2.2  Vortex  ring  impinging  on  a  wall 

To  demonstrate  the  ability  of  the  proposed  AMR  formulation  to  accurately 
represent  vorticity  dynamics,  the  simple  yet  very  challenging  problem  of  a  vortex  ring 
impinging  on  a  wall,  is  considered.  Due  to  the  importance  of  wall-vortex  interactions 
in  many  technological  applications  there  is  variety  of  reference  data  in  the  literature 
[104,  28,  73,  96].  In  the  present  work,  the  author  will  consider  a  setup  analogous  to 
the  one  reported  in  [96] ,  where  a  vortex  ring  is  generated  in  the  center  of  the  x  —  y 
plane,  and  at  a  distance  Zo  =  I.SDq  from  the  wall  as  shown  in  Fig.  3.9.  The  vortex 
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Figure  3.9:  Computational  setup  for  the  vortex  ring  impinging  on  a  wall 
ring  was  generated  by  introducing  an  impulsive  body  force  of  the  form: 


-AoT{t)F{z)H{r), 

T(t)  = 

0.5[1  -|-  tanh  Q;(r  —  t 

F(z)  = 

0.5[1  -|-  tanh/3(R^  — 

//(r)  = 

0.5[1  -|-  tanh7(C'r  — 

to\  and  z'  = 

\z  —  Zoj-  Setting  Ag 

(3.16) 


r  =  0.04,  {3  =  100,  =  0.1,  7  =  100  and  Cr  =  0.5,  the  resulting  vortex  ring 

had  Re  =  UoDo/p  =  570,  where  Uo  and  Do  are  its  initial  diameter  and  self-induced 
translation  velocity  and  u  is  the  kinematic  viscosity  of  the  fluid.  This  is  lower  than 
the  value  of  Rco  =  645,  reported  in  [96].  A  closer  match  of  the  initial  conditions  was 
not  possible  due  to  the  fact  that  some  of  their  forcing  parameters  were  not  listed. 

The  AMR  grid  was  adaptively  refined  in  areas  of  high  velocity  gradients.  Five 
levels  of  refinement  were  used  and  the  total  number  of  points  was  of  the  order  of 
2  X  10®.  The  grid  adaptation  was  done  every  10  timesteps  using  the  modulus  of 
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the  vorticity  field,  |a;|,  as  a  test  variable.  The  grid  was  refined  in  leaf  blocks  where 
|a;|  >  5.5Uo/Lo  anywhere  within  the  block,  and  derefined  if  <  A.OUo/ Lg.  Periodic 
boundary  conditions  are  applied  in  the  x  and  y  directions,  and  non  slip  at  the  top 
and  bottom  walls  in  the  direction  (see  Fig.  3.9).  The  AMR  solution  is  compared 
to  a  computation  of  exactly  the  same  problem  using  a  single-block,  finite-difference 
solver  [9].  The  single-block  grid  in  the  latter  case  utilizes  256  x  256  x  128  grid  cells, 
and  is  uniform  in  the  x  and  y  directions,  while  it  is  stretched  in  the  wall-normal 
direction,  ;2.  The  resulting  cell  size  in  the  near  wall  area  was  comparable  to  the  one 
at  the  highest  refinement  level  in  the  AMR  computation. 

In  Fig.  3.10,  vorticity  contours  are  shown  for  both  the  single-block  and  AMR 
calculations.  The  vorticity  normal  to  the  y  —  z  plane  is  shown  at  four  different  times 
during  the  calculation.  In  the  AMR  case  the  block  boundaries  are  also  shown  in 
the  figure.  The  similarity  between  the  two  computations  is  striking.  As  the  primary 
vortex  approaches  the  wall  there  is  an  increase  of  its  radius,  which  is  accompanied 
by  a  noticeable  decrease  in  the  core  size  (i.e.  Figs.  3.10a-c).  The  boundary  layer 
generated  underneath  the  primary  vortex  thickens  due  to  the  adverse  pressure  gra¬ 
dient  in  the  radial  direction  (Fig.  3.10b-f),  and  eventually  the  accumulated  vorticity 
pinches  off  forming  a  secondary  vortex  (Fig.  3.10c-g). 


57 


0.( 


IHHirK;; 


IIHJMfUNjWSBH 


Mwmnm 

OE^ilLZIUl 


Secondary 


IwtiiaBIIiMI 

lElSHOnil 


Figure  3.10:  Vorticity  isolines  at  an  y  —  z  plane  for  the  case  of  the  vortex  ring 
impinging  to  a  wall.  Forty  hve  ui^Ro/Uo  contours  from  —50  to  50  are  used.  Left 
side  is  from  the  single-block  calculation, and  the  right  side  from  the  AMR.  (a),(e) 
t  =  1.3;  (b),(f)  t  =  1.5;  (c),(g)  t  =  2.1;  (d),  (h)  t  =  2.5. 


(b) 


Figure  3.11:  Isosurfaces  of  Q  for  the  case  of  a  vortex  impinging  to  a  wall  at  three 
different  times.  The  AMR  block  distribution  is  shown  at  a.n  y  —  z  plane,  (a) 
t  =  2.0;  (b)  t  =  2.7;  (c)  t  =  4.4.  (i),  (ii),  and  (Hi)  indicate  the  primary,  secondary 
and  tertiary  vortices. 

A  three-dimensional  view  of  primary  and  secondary  vortex  structures  can  be 
seen  in  Figure  3.11,  where  isosurfaces  of  the  second  invariant  of  the  velocity  gradient 
tensor,  Q  =  —l/2{dui/dxjduj/dxi),  are  plotted.  The  secondary  vortex  orbits  the 
primary  vortex  (Fig.  3.11a),  and  later  undergoes  a  ’buckling’  instability  as  it  is 
compressed  by  the  primary  vortex  (Fig.  3.11b).  The  secondary  vortex  breaks  into 
smaller  structures  that  are  advected  under  the  primary  vortex  and  affect  the  vorticity 
generation  near  the  wall.  This  results  in  small  slender  structures  dominated  by  radial 
vorticity  (see  Fig.  3.11c).  At  the  same  time  a  third  vortex  pinches-off  and  also  starts 
orbiting  the  primary  one.  This  structure  at  this  Rco  is  fairly  stable,  and  stays  in 
the  proximity  of  the  primary  vortex  for  the  rest  of  the  calculation  (Fig.  3.11c). 

In  Figure  3.12  the  trajectories  of  the  primary  and  secondary  vortex  centers 
are  shown  for  both  computations.  The  numerical  results  from  Swearingen  et  al.  [96] 
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Figure  3.12:  Trajectories  of  the  centers  of  the  primary  and  secondary  vortices  at 
Rco  =  570.  Symbols  are  from  the  uniform  grid  and  lines  from  the  AMR  grid 
computation.  •  primary  vortex;  o  secondary  vortex;  green  line  is  the  trajectory  of 
the  primary  vortex  from  [96]  at  Rcq  =  645. 

for  a  similar  setup  at  Rco  =  645  are  also  included  for  comparison.  The  agreement 
between  the  single-block  and  AMR  calculations  is  excellent  and  the  trajectories 
of  both  vortices  are  identical  indicating  that  the  proposed  AMR  strategy  properly 
captures  the  vorticity  dynamics.  The  agreement  with  the  results  reported  in  [96]  is 
also  good.  Small  discrepancies  near  the  wall  are  primarily  due  to  the  ambiguity  in 
the  dehnition  of  the  vortex  center  as  well  as  the  differences  in  the  initial  conditions. 
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3.2.3  Fluid- Structure  interaction  of  two  falling  plates 

To  test  the  accuracy  of  the  approach  in  fluid-structure  interaction  problems 
the  case  of  two  falling  plates  is  considered  (see  Fig.  3.13).  In  this  case,  Eq.  (2.3) 
governing  the  dynamics  of  each  plate  can  be  reduced  to  the  following  form: 
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with  =  [x(t)  y(t)  tl2  =  di-  ^(t))  y(t)  are  the  coordinates  of  the  center  of 

mass  of  the  plate  in  the  x  and  y  directions  respectively,  and  9{t)  its  orientation  angle 
(see  Fig.  3.13a).  fx,  fy,  are  the  corresponding  hydrodynamics  forces  acting  on  the 
plate,  and  M^,  is  their  moment  with  respect  to  the  center  of  mass.  Eq.  (3.17)  has 
been  made  dimensionless  using  the  chord  length,  c,  the  mean  descent  velocity  of  plate 
2,  Ur,  and  the  fluid  density,  pf,  as  reference  variables.  The  thickness  of  each  plate 
is  10%  of  the  chord  length,  and  their  density  is,  ps  =  5.1p/.  Also  g  =  3.0,  m  =  0.5, 
and  lo  =  4.2  x  10“^.  The  resulting  Reynolds  number  is  Re  =  Urc/v  =  200.  The 
computational  domain  together  with  the  initial  conditions  are  shown  in  Fig.  3.13b. 
Both  plates  are  initially  at  rest.  To  evaluate  the  accuracy  of  the  proposed  AMR 
scheme  two  separate  simulations  are  considered.  First,  a  single  level  calculation 
using  a  uniform  grid  of  1024  x  1280  nodes  with  spacing  Ax  =  Ap  =  0.0078Lr  is 
conducted.  Then,  an  AMR  computation  with  four  levels  of  refinement  is  carried  out. 
Mesh  adaptivity  in  this  case  was  guided  by  two  criteria:  i)  the  presence  or  not  of  a 
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Figure  3.13:  (a)  Variables  describing  two-dimensional  rigid  body  motion  for  the 
falling  plates,  (b)  Domain,  boundary  and  initial  conditions  for  the  two  falling  plates 
problem. 

rigid  body  in  a  grid  block,  ii)  the  vorticity  modulus,  |u;|,  as  in  the  example  given  in 
the  previous  section.  As  a  result  the  highest  rehnement  level  always  surrounds  the 
plates,  as  well  as  the  areas  of  high  velocity  gradients  in  their  wake.  Note  that  grid 
cell  size  at  the  highest  rehnement  level  is  the  same  as  the  one  used  in  the  uniform 
grid  calculation.  Both  computations  were  advanced  in  time  for  22  x  10^  steps,  with 
At  =  2.0  X  10-h 

The  trajectories  of  the  centers  of  mass  for  both  bodies  as  a  function  of  time, 
and  phase  diagrams  of  positions  and  velocities  for  the  horizontal  and  vertical  degrees- 
of-freedom  are  shown  in  Figs.  3.14a,b  respectively.  The  agreement  with  the  single 
block  computation  is  excellent  demonstrating  the  hdelity  of  the  proposed  AMR 
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(a)  (b) 

Figure  3.14:  (a)  Position  of  center  of  mass  for  each  plate  as  a  function  of  time,  (b) 
Phase  diagrams  for  xi^t)-  Symbols  are  from  the  AMR  computation  and  lines  from 
the  uniform  grid  computation.  •  plate  1;  +  plate  2.  Also,  resulting  trajectories  and 
phase  diagrams  for  plate  2  with  coarser  resolution  uniform  grids  are  plotted:  (-  -) 
512  X  640  grid,  (-.-)  256  x  320  grid. 

formulation.  The  trajectories  and  phase  diagram  obtained  for  plate  2  with  two 
uniform  grids  of  coarser  resolution  than  the  one  obtained  with  4  rehnement  levels 
are  also  plotted  for  comparison.  Grid  convergence  of  the  solution  is  observed  as 
the  uniform-grid  resolution  is  increased.  We  see  here  the  effectiveness  of  the  AMR 
formulation  in  representing  the  uniform,  fine  grid  solution,  which  is  not  achievable 
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using  coarser  uniform  meshes. 

A  frame-by-frame  comparison  of  the  vorticity  held  in  the  two  computations  is 
also  shown  in  Fig.  3.15  with  very  good  agreement  as  well.  It  should  be  mentioned 
that  the  proposed  FSI  treatment  has  been  tested  on  a  single  grid-block  setting  by 
Yang  et  ah  [114],  on  different  incompressible  how  problems,  and  found  second  order 
accurate  in  space  and  time.  It  is  interesting  to  note  the  transition  from  steady 
buttering  fall  to  tumbling  for  plate  2  (see  Fig.  3.15d,f).  Tumbling  motion  has  been 
observed  by  Andersen  et  al.  [4]  for  a  single  falling  plate  under  similar  how  conditions, 
and  usually  occurs  for  Re  >  100  and  Jo  >  3  x  10“^  [18].  For  plate  1  transition  from 
a  buttering  fall  to  tumbling  is  limited  by  the  wake  of  plate  2. 

3.2.4  Three  dimensional  example:  Sphere-wall  collision 

To  investigate  the  robustness  and  accuracy  of  the  proposed  method  in  three 
dimensional  conhgurations  we  performed  computations  of  a  rigid  sphere  bouncing 
oh  a  wall.  Problems  involving  collisions  between  immersed  bodies  are  particularly 
challenging  for  direct-forcing  schemes,  since  the  presence  of  two  or  more  Lagrangian 
markers  from  diherent  bodies  in  the  proximity  of  the  same  Eulerian  grid  cell  is 
usually  the  source  of  ambiguity.  The  proposed  forcing  scheme  treats  such  situations 
in  a  robust  manner  without  the  need  for  special  treatments. 

The  particular  conhguration  we  selected  has  applications  to  particulate  hows, 
and  a  number  of  experimental  (i.e.  [31],  [48])  and  numerical  (i.e.  [6])  results  are 
available  in  the  literature  for  comparison.  The  dominant  parameter  in  the  collision 
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Figure  3.15:  Comparison  between  uniform  (left  side)  and  AMR  (right  side)  grid 
computations  for  the  case  of  the  falling  plates.  Vorticity  isolines  are  shown  at  (a)  & 
(d)  t  =  4.6;  (b)  &  (e)  f  =  5.2;  (c)  &  (f)  t  =  5.8. 
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process  is  the  Stokes  number,  St  =  l/9{pb/pf)Re,  where  pb  and  pf  are  the  particle 
and  fluid  densities  respectively,  and  Re  is  the  Reynolds  number  based  on  the  particle 
diameter,  D,  and  the  translational  velocity,  Uf,  an  instant  before  impact.  For  low 
values  of  the  Stokes  number  [St  <  10)  no  rebound  will  occur,  even  if  the  dry 
restitution  coefficient,  Cdry,  is  different  from  zero.  For  St  >  10  rebound  occurs,  and 
the  total  restitution  coefficient,  er,  is  lower  than  edry  For  large  values  of  the  Stokes 
number  {St  >  500)  the  total  restitution  coefficient,  e-r,  approaches  edry  If  ^dry  =  0, 
then  the  Reynolds  number  is  the  only  dominant  parameter. 

Below  we  will  present  results  from  two  different  configurations:  i)  a  case  where 
no  rebound  occurs,  which  resembles  the  conditions  in  the  experiments  by  Fames  and 
Dalziel  [31];  ii)  a  case  where  rebound  is  allowed,  which  resembles  the  conditions  in 
the  axisymmetric  Navier-Stokes  computations  by  Ardekani  et.  al  [6].  In  figure  3.16a 
a  sketch  of  the  computational  domain  is  shown.  Both  the  sphere  and  the  wall  are 
immersed  into  a  locally  refined  Cartesian  mesh  [11],  and  are  represented  by  an 
unstructured  Lagrangian  grid  with  2.7  x  10^  and  2.9  x  10^  markers  respectively  (see 
figure  3.16b).  The  Eulerian  grid  is  arranged  in  a  way  that  the  resolution  around  the 
sphere  is  Ax  =  O.OIH.  In  both  cases  the  horizontal  displacements  and  rotations  are 
constrained  and  the  vertical  displacement,  Zs{t),  is  governed  by: 

msZs{t)  =  -msQ  +  fz{zs,  Zs,  t),  (3.18) 

where  m*  is  the  mass  of  the  sphere,  g  is  the  acceleration  of  gravity  and  fz  is  the  hydro- 
dynamic  force  on  the  sphere  in  the  vertical  direction.  The  Navier-Stokes  equations 
governing  the  dynamics  of  the  fluid,  and  equation  (3.18)  governing  the  dynamics 
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Figure  3.16:  Sphere- wall  interaction  problem,  (a)  computational  setup;  (b)  La- 
grangian  and  Eulerian  grids. 

of  the  sphere  are  solved  as  a  coupled  system  using  the  predictor-corrector  strategy 
proposed  in  [114]. 

Initially  the  sphere  is  located  at  a  distance  of  4.3/1  from  the  horizontal  wall, 
and  is  impulsively  started  to  reach  a  velocity  corresponding  to  an  initial  Reynolds 
number  of  Rci  =  510.  The  density  ratio  is  hxed  to  Ps/Pf  =  3.2.  The  resulting 
Reynolds  number  just  before  impact  is  Re  =  830  and  the  corresponding  Stokes 
number  is  St  =  295.  Contact  is  assumed  to  take  place  when  sphere  and  floor  are 
within  a  distance  of  one  cell  size.  In  the  experiments  by  Eames  and  Dalziel  [31]  the 
Reynolds  number  before  impact  was  Re  =  850,  but  the  motion  of  the  sphere  was 
prescribed  and  no  bounce  was  allowed  to  occur.  To  simulate  these  conditions  in 
our  computations  the  dry  restitution  coefficient,  edry,  was  set  to  zero  (see  equation 
(3.19)  below). 
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In  Figure  3.17  dye  flow  visualizations  form  the  experiments  in  [31]  are  directly 
compared  to  our  computations,  where  the  flow  patterns  are  visualized  using  az¬ 
imuthal  vorticity  isolines.  The  computed  flow  patterns  are  in  very  good  qualitative 
agreement  with  the  experiment.  Direct  quantitative  comparisons  are  not  possible 
due  to  the  fact  that  vorticity  and  scalars,  such  as  dye,  do  not  have  the  same  dynam¬ 
ics  as  a  result  of  their  different  diffusivities  and  the  absence  of  vortex  stretching  in 
the  case  of  scalars.  As  the  sphere  approaches  the  wall  the  detached  shear  layers  and 
the  small  recirculating  areas  behind  the  sphere  are  evident  in  both  experiments  and 
computations  (figure  3.17a,  6).  Just  after  impact  the  vorticity  in  the  shear  layers 
moves  towards  the  wall  generating  a  layer  of  vorticity  with  opposite  sign  on  the 
sphere’s  surface  (see  figure  3.17c).  As  soon  as  this  layer  separates  (figure  3.17d), 
a  vortex  dipole  is  formed  and  moves  away  from  the  sphere  (figures  3.17e,/).  By 
locating  the  center  of  these  vortices  in  the  computations  and  the  experiments  we 
found  that  in  their  trajectories  through  time  are  always  within  5%. 

Next,  we  considered  the  bouncing  sphere  problem.  The  dry  restitution  coef¬ 
ficient  was  the  one  used  in  the  axisymmetric  calculations  of  Ardekani  et.  al  [6], 
Gdry  =  0.97,  which  is  typical  for  steel-sphere  and  glass-wall  collisions.  The  collision 
process  starts  when  the  distance  between  the  particle  and  the  wall  is  equal  to  the 
roughness  height,  hr-  We  assume  that  rough  surface  has  a  negligible  effect  on  the 
viscous  force  until  the  gap  between  the  smooth  portions  of  surfaces  becomes  equal 
to  the  size  of  largest  roughness  element,  hr-  This  is  also  the  moment  the  impact  is 
assumed  to  occur.  Details  can  be  found  in  [6].  Just  after  the  collision  we  define  a 


Figure  3.17:  Sphere-wall  interaction  with  Cdry  =  0.0.  The  left  half  in  all  hgures  are 
azimuthal  vorticity  isolines  from  the  present  computations,  and  the  right  half  are 
snapshots  from  the  dye  visualizations  in  [31]  at  Re  =  850.  (a)  — r/,  (b)  0,(c)  T/,(d) 
2r/,  (e)  3r/,  and  (f)  4r/,  where  Tf  =  D/Uf. 

new  set  of  initial  conditions  for  equation  (3.18)  as  follows: 


^s2  ^s2  ^dry^sl  (3.19) 

where  Zgi,  Zgi  and  Zs2,  Zs2  are  the  sphere’s  vertical  position  and  velocity  before  and 
after  the  impact  respectively. 

In  numerical  simulations  of  contact  problems  it  is  important  that  the  lubri¬ 
cation  layer  between  the  bodies  is  resolved.  In  all  our  computations  the  surface 
roughness,  which  practically  determines  how  close  the  bodies  can  come,  and  numer¬ 
ical  resolution  were  selected  in  a  way  that  a  minimum  of  5  —  6  Eulerian  grid  points 
were  present  between  the  bodies  during  impact.  We  hrst  considered  the  same  case 
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as  in  the  above  example  {Re  =  830  and  St  =  295)  with  a  non-zero  restitution  co¬ 
efficient,  Cfiry  =  0.97,  in  order  to  compare  the  vortex  dynamics  with  the  no-bounce 
case.  In  hgure  3.18  azimuthal  vorticity  isolines  of  the  no-bounce  (left)  and  bouncing 
sphere  (right)  are  shown.  Just  before  impact  (hgure  3.18a)  the  how  for  both  cases 
is  identical  since  we  start  from  the  same  initial  conditions.  At  a  later  time  and  after 
the  hrst  impact  (3.186,  c),  the  layers  of  vorticity  with  alternating  sign  that  were  ob¬ 
served  in  the  case  with  Cdry  =  0.0  can  also  be  seen  in  the  bouncing  sphere  problem. 
In  the  latter  case,  however,  the  primary  vortex  originating  in  the  wake  is  weaker 
and  the  upward  motion  of  the  sphere  causes  the  shear  layer  at  the  surface  to  roll-up 
into  a  strong  secondary  vortex.  As  the  downward  motion  of  the  sphere  starts  the 
secondary  vortex  pinches-oh  (see  hgure  3.18d)  and  by  the  time  the  second  bounce 
occurs  it  is  dissipated.  As  a  result  the  wake  and  secondary  vortices  do  not  form  the 
dipole  structure  seen  in  the  no  bounce  case. 

Ardekani  and  Rangel  [6]  dehned  a  total  restitution  coefficient,  ct  =  Ua/Uf, 
where  Ua  is  the  velocity  of  the  sphere  at  tUf/D  =  0.07  after  the  impact  time,  tc, 
which  measures  the  dissipative  ehect  of  the  huid,  as  it  is  drained  and  subsequently 
reenters  the  layer  between  sphere  and  wall.  For  the  present  case  we  found  ct  =  0.63. 
A  direct  comparison  with  the  computations  reported  in  [6],  where  they  reported 
ct  =  0.92,  is  not  possible  because  of  the  diherences  in  the  Reynolds  numbers.  Our 
ct  value,  however,  is  consistent  with  the  trend  reported  in  [6]  where  a  decrease  in 
ct  was  observed  with  increasing  Reynolds  number  and  constant  Stokes  number.  In 
particular  they  found  a  decrease  of  5.0%  for  er  at  St  =  301  when  Re  increases 
from  35  to  162.  To  further  verify  the  accuracy  of  our  formulation  we  also  conducted 
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a  computation  that  closely  matches  the  low  Reynolds  number  conditions  in  the 
simulation  by  Ardekani  and  Rangel  [6].  The  Reynolds  Number  before  impact  was 
Re  =  76.8  and  the  Stokes  number  St  =  299.  In  this  case  we  computed  an  ct  =  0.91, 
which  is  in  very  good  agreement  with  the  reference  results  of  ct  =  0.92. 


(c) 


(f) 

Figure  3.18:  A  comparison  of  Cdry  =  0.00  (left)  and  Cdry  =  0.97  (right)  computations 
for  the  sphere-wall  interaction  at  Re  =  830.  (a)  — O.lr/,  (b)r/,  (c)2r/,  (d)  3r/,  (e) 
3.5r/,  and  (f)  4.1r/,  where  Tf  =  D/Uf. 


T  i 
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Chapter  4 


Large  eddy  simulation  for  discontinuous  grids 

In  the  introductory  Chapter  we  described  a  set  of  features  of  LES  which  make 
its  combination  with  S-AMR  methods  non  trivial.  The  grid  discontinuities  intro¬ 
duced  in  AMR  methods  may  lead  to  numerical  errors,  in  some  cases  signihcant  ones. 
This  is  especially  true  if  non-dissipative,  centered  spatial  discretization  schemes  like 
the  one  previously  proposed  are  used.  In  the  following  section  we  give  a  formal 
introduction  to  eddy  viscosity  models  in  large  eddy  simulations.  We  then  develop  a 
simple  strategy  to  vary  the  hlter  size  for  hltered  variables  around  grid  discontinu¬ 
ities.  In  order  to  evaluate  the  sources  of  error  and  develop  a  strategy  to  minimize 
these,  we  will  study  in  Section  4.2  the  flow  of  spatially  decaying  isotropic  turbulence 
advected  across  a  rehnement  interface.  In  this  study  we  see  that  using  LES  together 
with  explicit  hltering  of  the  advective  terms  [60]  leads  to  an  effective  control  of  alias¬ 
ing  and  interpolation  errors  when  turbulence  travels  across  a  derehnement  interface. 
This,  strategy  is  used  in  Section  4.3  in  the  simulation  of  transitional  flow  around  a 
sphere  at  Re  =  10000. 

4.1  Mathematical  formulation 

We  employ  overbars  to  denote  the  resolved  quantities  on  LES  computations. 
We  consider  the  case  where  the  SGS  stresses,  Tij  in  the  momentum  equations  (2.1), 
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are  parameterized  using  standard  eddy-viscosity  models: 

2  — 

'^ij  ij  1  (^'1) 

where  ut  is  the  eddy  viscosity  dehned  as: 


ut  =  CA}\S\; 


(4.2) 


C*  is  a  model  coefficient,  A/  the  hlter  width,  and  Sij  the  resolved  strain-rate  tensor: 


C  f  ^  _L  ^ 
2  \dxi  dxj 


(4.3) 


In  this  work  we  nse  two  different  models  to  compnte  vt'-  the  standard  Smagorinsky 
model  [92,  56],  and  the  Lagrangian  dynamic  eddy- viscosity  (LDEV)  model  proposed 
by  Menevean  et  al.  [66].  In  the  former  the  model  coefficient  is  a  constant  set  to 
C  =  0.65^,  while  in  the  latter  it  is  dynamically  compnted  dnring  the  calculation: 


where 


1  {CjjMjj) 

2  ’ 


Mji  =  A  f 


•■ij 


S 


Aij  XliUj  Uj^Uj^ 


(4.4) 


(4.5) 

(4.6) 


Note  that  7  is  the  test-filter  operator,  representing  a  filter  with  width  At  larger 
than  the  grid  filter  A/,  and  a  =  At/ A/  is  the  ratio  between  the  test  and  grid 
filter-widths.  The  numerator  and  denominator  in  equation  (4.4),  are  averaged  over 
particle  pathlines.  Details  on  the  present  implementation  of  the  Lagrangian  model 
can  be  fonnd  in  [90]. 
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Any  filtered  quantity,  /  and  /,  is  obtained  by  recursive  application  of  a  one- 
dimensional  filter  in  the  three  coordinate  directions;  the  one-dimensional  hlter  is 
dehned  as: 


j+M 

f^=  Y,  (4.7) 

m=j—M 


(and  similarly  for  /).  The  hlter-width  corresponding  to  above  operation  is  given  by 

|59]: 

/  M  \ 

A/  =  A  12  E 

\  j=-M  ) 

where  A  is  the  grid  size.  In  our  computations  we  set  =  a/4/3  and  a  = 


At/Af  =  Vq. 


4.1.1  Explicit  filtering  of  the  non-linear  term 

When  explicit  hltering  of  the  non-linear  term  is  used  in  LES,  the  hltered 
velocity  product  in  the  advective  terms  of  the  momentum  equations  is  decomposed 
as  [54]: 

UiUj  =  UiUj  -|-  {uiUj  —  UiUj)  (4.9) 

where  rh  =  UjUj  —  UiUj  is  the  new  subgrid  stress  tensor  to  be  modeled.  The  modeling 
strategy  used  here  is  the  same  as  discussed  in  the  previous  Section,  and  the  same 
hlters  are  applied.  The  secondary  filtering  operation  has  the  effect  of  inhibiting  the 
generation  of  frequencies  higher  than  the  grids  characteristic  wavenumber  [103],  thus 
reducing  the  numerical  error  at  the  higher  wave-numbers  included  by  the  mesh  (this 
is  particularly  important  in  the  framework  of  the  LDEV  model,  as  the  parameter 
C  and  the  SGS  stresses,  depend  intrinsically  on  the  accuracy  at  which  the  smallest 
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resolved  scales  are  computed).  However,  explicit  filtering  as  defined  by  equation 
(4.9)  is  not  in  general  Galilean  invariant  with  respect  to  translation.  Lund  [60] 
showed  that  the  invariance  error  is  reduced  if  the  hlter  used  becomes  closer  to  the 
sharp  Fourier  cutoff;  adding  a  scale-similarity  type  term  to  the  model  also  leads  to 
a  Galilean-invariant  formulation  [81].  The  present  study  is  limited  to  eddy- viscosity 
models  with  trapezoidal  Liters  and  Galilean  invariance  is,  therefore,  not  enforced. 
The  error,  however,  is  proportional  to  the  to  velocity  of  the  moving  reference  frame 
[60],  which  in  our  case  is  zero. 

4.2  Spatially  decaying  isotropic  turbulence  past  a  refinement  inter¬ 
face 

To  quantify  the  numerical  errors  arising  on  LES  in  the  vicinity  of  grid  refine¬ 
ment  interfaces,  and  explore  ways  of  reducing  them,  we  study  a  simple  flow  con¬ 
figuration,  without  some  of  the  complexities  present  in  wall-bounded  flows  (mean 
shear,  turbulence  anisotropy,  etc.).  We  examine  homogeneous  isotropic  turbulence 
advected  by  a  uniform  velocity,  Gc,  across  a  fine-coarse  grid  interface,  concentrating 
on  the  case  in  which  the  interface  is  normal  to  the  main  advection  direction.  The 
development  of  the  turbulence  in  the  vicinity  and  downstream  of  the  interface  are 
examined  to  evaluate  the  errors  introduced  by  the  grid  discontinuity,  and  the  dis¬ 
tance  required  for  a  return  to  equilibrium.  We  also  investigate  the  effect  of  explicit 
hltering  the  advective  term  in  the  momentum  equations  as  a  means  of  separating  the 
hlter  size  from  the  hnite-difference  grid  size,  and  to  control  the  frequency  content 
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Case 

Grid 

Filter  Variation 

C 

48^  +  48^ 

N/A 

F 

96^  -h  96^ 

N/A 

FC-S 

96^  +  48^ 

Sharp 

CF-S 

48^  +  96^ 

Sharp 

FC-V 

96^  +  48^ 

Linear  variation  over  l.QLr 

Table  4.1:  Summary  of  Lagrangian  Dynamic  inflow-outflow  computations, 
and  numerical  errors  of  the  solution  close  to  the  interface. 

4.2.1  Computational  setup 

We  conducted  simulations  of  spatially  decaying  homogeneous  isotropic  tur¬ 
bulence.  A  typical  conhguration  is  shown  in  Figure  4.1  (left).  The  computational 
domain  consists  of  two  adjoining  blocks,  each  with  dimensions  27r  x  27r  x  27r  in  the  x, 
y  and  directions  respectively.  Each  block  is  resolved  by  24^,  48^  or  96^  grid  points; 
the  mesh  is  uniform  and  the  grid  size  is  the  same  in  all  directions.  The  two  blocks 
have  the  same  resolution  for  the  single-grid  calculations,  while  in  the  two-level  cases 
the  grid  size  is  discontinuous,  with  a  factor  of  two  ratio  in  each  direction  (he.,  we 
go  from  a  48^  grid  to  a  96^  one  in  the  coarse-to-fine  cases,  from  a  96^  grid  to  a  48^ 
one  in  the  £ne-to-coarse  ones).  The  refinement  jump  (the  coarse-fine  grid  interface) 
is  located  at  the  center  of  the  domain  {x  =  0),  normal  to  the  x  direction.  This 
interface  defines  two  sets  of  calculations,  coarse-to-fine  (CF)  as  in  Figure  4.1,  or 
fine-to- coarse  (FC).  Tables  4. 1-4. 2  show  the  sets  of  simulations  carried  out  for  the 
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Case 

Grid 

Filter  Variation 

C 

48^  +  48^ 

N/A 

F 

96^  +  96^ 

N/A 

FC-S 

96^  +  48^ 

Sharp 

CF-S 

48^  +  96^ 

Sharp 

FC-V 

96^  +  48^ 

Linear  variation  over  l.GL^ 

CF-V 

48^  +  96^ 

Linear  variation  over  l.GL^ 

Table  4.2:  Summary  of  Smagorinsky  inflow-outflow  computations. 


Outflow 


Figure  4.1:  Setup  of  the  numerical  experiments. 

LDEV  and  Smagorinsky  models  respectively. 

To  generate  the  inflow  conditions,  computations  of  forced  isotropic  turbulence 
at  the  same  resolution  as  the  block  adjacent  to  the  inflow  plane  were  hrst  conducted. 
The  linear  turbulence  forcing  scheme  proposed  by  Lundgren  [61]  and  further  studied 
by  Rosales  and  Meneveau  [85]  was  applied,  using  a  value  of  mean  dissipation  per 
unit  mass,  £  =  ^AUf/Lr,  and  a  Reynolds  number  Rcr  =  10,000.  Planes  from  this 
simulation  were  stored  on  disk  and  used  as  inflow  conditions  in  the  inflow-outflow 
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Case 

Grid 

Filter  Variation 

C 

243  +  2  X  243 

N/A 

F 

48^  +  48^ 

N/A 

FC-S 

48^  +  2  X  243 

Sharp 

CF-S 

243  +  2  X  48^ 

Sharp 

FC-V 

48^  +  2  X  243 

Linear  variation  over  1.6L,. 

CF-V 

243  +  2  X  48^ 

Linear  variation  over  l.GL^ 

Table  4.3:  Summary  of  Lagrangian  Dynamic  inflow-outflow  computations  on  lower- 
resolution  grids. 

computations  (see  Figure  4.1).  In  order  to  maintain  consistency,  the  same  SGS 
model  was  applied  for  the  periodic  and  inflow-outflow  simulations.  Note,  however, 
that  the  results  of  the  forced  simulations  were  found  to  be  quite  sensitive  to  the 
model  and  grid  resolution;  thus,  at  the  inflow  we  had  slightly  different  values  of 
turbulent  kinetic  energy,  depending  on  the  model  and  resolution  used.  A  specified 
convective  velocity,  Uc  =  0.4114-,  is  used  to  advect  turbulence  in  x  direction.  This 
value  was  chosen  based  on  two  considerations:  hrst,  we  wanted  a  ratio  Uc/urms 
representative  of  the  conditions  encountered  in  a  wall-bounded  flow,  especially  near 
the  wall  where  the  turbulent  eddies  are  smaller;  this  would  result  in  Uc/urms  —  4. 
On  the  other  hand,  we  wanted  to  use  as  short  a  domain  as  possible,  to  reduce  the 
computational  cost  of  the  simulation.  As  a  result,  we  chose  Uc/urms  —  2.5.  A 
convective  condition  [74]  is  used  at  the  outflow  and  periodic  conditions  in  all  other 
directions. 
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x/L^ 


x/L^ 


Figure  4.2:  Filter  width  variation. 

For  each  model  we  performed  single-grid  calculations  on  coarse  and  fine  grids, 
as  well  as  cases  in  which  the  grid  and  filter  are  changed  abruptly  by  a  factor  of  two 
in  each  direction  (FC-S  and  CF-S),  and  cases  in  which  the  filter  width  was  linearly 
increased  or  decreased  (on  the  fine-grid  side)  from  the  value  corresponding  to  the 
hne  mesh  to  that  corresponding  to  the  coarse  one  (CF-V  and  FC-V  cases).  The 
variation  took  place  over  a  distance  l.bL^  (comparable  to  the  integral  scale  of  the 
flow);  variations  over  shorter  distances  did  not  prove  to  be  effective.  Figure  4.2 
illustrates  the  variation  of  A f  (normalized  by  the  grid  size  on  the  fine- resolution 
side)  for  the  various  simulations. 

Note  that  the  use  of  a  filter-width  different  from  the  grid  size  is  ambiguous, 
unless  the  non-linear  term  is  hltered  explicitly.  When  the  Smagorinsky  model  is 
used,  for  instance,  it  is  not  possible  to  distinguish  a  change  in  Aj  from  a  change 
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in  C.  With  dynamic  models,  moreover,  A/  does  not  enter  the  equations  directly: 
what  matters  are  the  ratio  a  =  At/Af,  and  the  weights  used  to  dehne  /.  We  take 
the  view  here  that  the  hlter-width  A f  can  be  assigned  a  priori,  independently  from 
the  grid  size.  The  coefficient  C  is  determined  from  the  flow  physics,  either  through 
integration  of  the  spectrum  for  the  Smagorinsky  model  (see  the  derivation  by  Lilly 
[56])  or  dynamically.  For  the  dynamic  model,  this  results  in  the  use  of  weighting 
coefficients  in  (4.7)  such  that  a  =  ^/6  independent  of  the  value  of  Aj  used. 

For  each  calculation,  data  at  several  locations  along  the  longitudinal  direction 
was  collected  in  the  course  of  more  than  120  computational  time  units.  Then  the 
statistics  were  computed  using  an  ensemble  of  no  less  than  100  time  snapshots,  and 
quantities  were  also  plane-averaged  for  each  x  location. 

4.2.2  Single-grid  calculations 

First,  we  performed  single-grid  calculations  with  a  coarse  and  a  hue  grid,  using 
the  Smagorinsky  and  LDEV  models.  The  coarse  grid  uses  48  points  in  each  direction 
for  each  block  (for  a  total  of  2  x  48^  =  2.2  x  10^  points)  while  the  hne  one  uses  96^ 
points  per  block,  for  a  total  of  1.8  million  points. 

In  these  calculations  the  mesh  is  cubic,  and  the  hlter-width  is  uniform;  commu¬ 
tation  errors  are,  therefore,  identically  zero.  As  expected,  slight  differences  between 
the  two  models  and  of  course  differences  between  the  calculations  with  the  coarse 
and  hne  grids  were  found.  The  explicit  hltering  of  the  advective  term  results  in  a 
larger  value  of  the  Taylor  microscale  A,  rehecting  the  fact  that  the  hltering  reduces 
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the  small-scale  energy  transfer  to  small  scales.  At  a;  =  0,  where  the  coarse-fine 
grid  interface  is  located  in  the  calculations  with  grid  refinement  or  coarsening,  the 
Reynolds  number  is  Re\  ~  800.  Also,  the  integral  length  scale  Ln  (obtained  by 
integrating  the  longitudinal  two-point  correlation)  is  approximately  O.SSL/j,  remain¬ 
ing  approximately  constant  along  the  streamwise  direction  and  between  the  different 
eddy-viscosity  models. 

4.2.3  Two-level  computations  with  the  LDEV  model 

First,  we  examine  the  results  obtained  with  the  LDEV  model.  Notice  that 
this  model  includes  memory  effects  through  the  integration  of  the  numerator  and 
denominator  of  (4.4)  along  Lagrangian  particle-paths.  Furthermore,  the  variable 
hlter-width  does  not  affect  the  eddy  viscosity  directly  {vt  depends  only  on  a),  but 
only  through  the  explicit  hltering  performed  at  the  test-hlter  scale. 

Figure  4.3  shows  isosurfaces  of  the  second  invariant  of  the  hltered  velocities 
gradient  tensor: 

^  ^  “2  ^  “2 

colored  by  the  value  of  the  eddy  viscosity  for  three  cases.  The  flow  behavior  in 
the  Coarse-to-Fine  (CF)  case.  Figure  4.3(a),  is  relatively  benign:  shortly  after  the 
interface  (denoted  by  the  dark  plane  at  a;  =  0)  we  observe  a  rapid  generation  of 
small  scales.  When  the  grid  is  suddenly  coarsened.  Figure  4.3(b),  on  the  other 
hand,  we  see  a  sudden  depletion  of  eddies  following  the  interface,  and  an  abrupt 
increase  in  their  length  scale.  When  the  hlter  is  smoothly  coarsened  on  the  hne-grid 
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side,  Figure  4.3  (^cj,  we  notice  an  increase  in  the  characteristic  eddy  length-scale  and 
an  increase  of  the  eddy  viscosity  even  upstream  of  the  interface,  and  a  smoother 
transition  to  the  coarse-grid  side. 

Figure  4.3  shows  that,  as  turbulence  is  advected  through  a  grid  discontinuity, 
the  flow  behavior  is  quite  different  depending  on  whether  the  grid  is  coarsened  or 
rehned.  The  changes  in  the  turbulence  structure  are  quite  dramatic,  and  one  should 
expect  them  being  reflected  in  the  flow  statistics. 

Figures  4.4  and  4.5  show  the  development  of  statistical  results.  It  should  be 
remarked  here  that  at  the  interface  between  hne  and  coarse  grids  the  interpolation 
is  used  to  exchange  information  between  hne  an  coarse  grids.  Furthermore,  because 
of  the  staggering,  the  control  points  where  the  velocity  components  are  dehned  on 
the  coarse  grid  do  not  coincide  with  those  on  the  hne  grid.  As  a  consequence,  while 
hrst-order  averaged  statistics  are  continuous  across  the  interface,  quantities  that 
are  quadratic  in  the  velocity  present  a  discontinuity.  Reducing  the  magnitude  of 
this  discontinuity  is  an  important  feature  of  a  successful  transition  methodology. 
For  the  CF  run  the  how  statistics  switch  between  the  coarse-  and  hne-grid  values 
reasonably  rapidly.  In  Figure  4.4,  the  longitudinal  spectra  at  various  x  locations  are 
shown.  These  are  the  transforms  of  the  two-point  correlation  of  the  in-plane  velocity 
components  with  displacement  in  the  direction  of  the  velocity  itself.  The  spectra 
immediately  after  the  coarse-hne  grid  interface  decay  too  steeply,  as  the  high  wave- 
numbers  initially  contain  no  energy.  They  are  hlled  rapidly  and  by  x/ Lr  =  0.7,  less 
than  one  integral  scale  downstream  of  the  interface,  the  spectrum  returns  to  the 
single-grid  value. 
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Figure  4.3:  Isosurfaces  of  Q,  colored  by  the  value  of  the  eddy  viscosity.  LDEV 
model,  (a)  CF-S;  (b)  FC-S;  (c)  FC-V.  The  flow  is  from  top  right  to  bottom  left. 
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Figure  4.4:  Longitudinal  spectra,  LDEV  model.  Each  spectrum  is  shifted  upward 
by  100  units  for  clarity. 

The  length  of  the  transition  downstream  of  the  grid  interface  depends  on  the 
quantity  examined,  but  is  generally  less  than  two  integral  scales  Ln  (recall  that 
Lii/Lr  —  0.85).  This  distance,  of  course,  is  strongly  dependent  on  Uc-  As  mentioned 
in  before,  the  value  used  here,  Uc  =  O.dlf/^,  corresponds,  at  the  interface,  to  a  mean 
velocity  that  is  about  2.5  times  the  rms  intensity  Urms]  in  boundary  layers  Uc/urms 
is  of  the  order  of  4  near  the  wall,  20  or  more  in  the  outer  layer;  these  values  would 
result  in  longer  transitions. 

Coarsening  the  grid  (EC  simulations)  introduces  a  much  stronger  perturbation. 
When  the  hlter  width  is  discontinuous  (FC-S  simulation)  we  observe  an  accumu¬ 
lation  of  energy  at  the  small  scales,  apparent  in  the  spectra  at  x  =  —0.25  and  0 
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Figure  4.5:  Streamwise  distribution  of  (a)  turbulent  kinetic  energy  K]  (b)  integral 


length  scale  Ln;  and  (c)  normalized  eddy  viscosity,  pt/v- 


LDEV  model. 
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(see  Figure  4.4).  The  energy  pileup  results  in  an  overshoot  of  the  turbulent  kinetic 
energy  K  =  (u'u')/2  immediately  before  the  interface  Figure  4.5(a).  We  also  ob¬ 
serve  discontinuities  in  the  integral  scale  and  the  eddy  viscosity  (normalized  by  the 
molecular  one). 

The  smoothly  increasing  hlter-width  (FC-V  simulation)  is  superior  to  the  dis¬ 
continuous  one:  since  the  LDEV  model  is  not  very  dissipative  (compared,  for  in¬ 
stance,  to  the  Smagorinsky  model  discussed  below),  the  increase  of  the  eddy  viscosity 
upstream  of  the  interface  is  needed  to  reduce  (but  not  completely  remove)  the  en¬ 
ergy  pileup  at  the  grid  interface  resulting  in  a  much  smoother  transition.  Even  in 
the  FC-S  case,  however,  despite  the  fact  that  the  energy  pileup  at  small  scale  is 
very  signihcant,  the  flow  adjusts  very  quickly  to  the  grid  coarsening,  and  returns  to 
the  single-grid,  uniform  hlter-width,  data  within  less  than  an  integral  scale  of  the 
interface. 

The  calculations  discussed  until  now  were  characterized  by  good  resolution  of 
the  integral  scales,  even  in  the  coarse-grid  size.  In  fact,  the  ratio  of  integral  scale 
L  =  je  (where  is  the  turbulent  kinetic  energy  and  e  the  dissipation)  to  the  grid 
size  is  approximately  15  on  the  48^  grid,  30  for  the  96^  one.  In  typical  wall-bounded 
how  calculations,  much  lower  resolution  is  generally  used:  in  the  buher  layer  and 
the  beginning  of  the  logarithmic  layer,  typically  5-10  grid  points  per  L  are  used, 
and  only  in  the  outer  layer  we  reach  resolutions  of  10-20  grid  points  per  integral 
scale.  To  represent  more  closely  the  resolution  levels  used  in  actual  calculations 
of  wall-bounded  hows,  we  performed  another  set  of  simulations  that  used  coarser 
meshes:  24  grid  points  were  used  in  all  directions  on  the  coarse  side  to  resolve  a 


271^  cube,  and  48  grid  points  on  the  fine  side  (see  Table  4.3).  Except  for  the  fine- 
level  calculation,  the  domain  after  the  interface  in  these  cases  was  extended  to  47r 
in  the  x  direction  to  reduce  the  effect  of  the  convective  outflow  boundary  on  the 
results.  In  the  24^  linearly  forced  simulations  used  to  obtain  inflow  data,  the  value 
of  dissipation  was  set  to  e  =  0.36U^ /Lr-  As  before,  in  all  two-level  calculations 
the  interface  was  located  at  a;  =  0,  and  the  convective  velocity  Uc  =  0.41f4.  The 
Reynolds  number,  based  on  A/,  at  the  interface  location  was  —  1500  for  the 
coarse  inflow-outflow  run. 

In  Figure  4.6  the  author  shows  the  average  non  dimensional  values  of  turbulent 
kinetic  energy,  longitudinal  integral  length  scale  and  normalized  turbulent  viscosity 
for  the  lower-resolution  simulations.  In  a  similar  trend  as  seen  in  the  past  section, 
the  Coarse-to-Fine  grid  and  filter- width  discontinuity  does  not  affect  the  statistics 
as  much  as  the  Fine-to-Coarse  one.  For  the  FC  cases,  again,  the  energy  pileup 
upstream  of  the  interface  is  corrected  by  smoothly  increasing  the  filter  width  more 
gradually.  The  sharp  FC  case  exhibits  a  larger  energy  pileup  zone  compared  to  the 
higher- resolution  simulation  of  Figure  4.5(a).  As  expected,  the  integral  length  scale 
is  Lii/Lr  =  1  a  value  higher  than  that  predicted  on  the  finer  meshes,  Ln/Lr  =  0.85. 
For  the  CF  case  the  varying  filter  in  the  fine-mesh  region  delays  the  return  to  the  fine 
grid  value  of  all  statistics:  the  increased  dissipation  associated  with  the  larger  filter- 
width  on  the  fine  grid  delays  the  generation  of  small  scale.  The  sharp  discontinuity 
CF-S  results  in  a  recovery  of  the  fine-grid  results  over  1.6Lr,  a  larger  distance  than 
observed  in  the  better  resolution  cases  discussed  above.  Similar  observation  may 
be  made  for  the  FC  cases:  the  smooth  filter-width  change  gives  improved  results. 
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and  the  length  required  for  recovery  from  the  grid  discontinuity  is  larger  than  in  the 
previous  case.  The  varying  hlter  CF  computation  exhibits  a  bump  in  the  curve, 
possibly  due  to  the  increase  of  [S'!  as  smaller  scales  are  generated,  while  the  hlter 
width  is  still  interpolated  from  the  coarse  side. 

4.2.4  Two-level  simulations  with  the  Smagorinsky  model 

Next  we  discuss  simulations  in  which  the  SGS  stresses  were  parameterized  us¬ 
ing  the  Smagorinsky  model.  The  Lagrangian  model  includes  memory  effects;  thus, 
the  eddy  viscosity  immediately  after  an  interface  is,  to  some  extent,  affected  by 
its  previous  history  and  recalls  the  fact  that  it  is  coming  from  a  region  of  coarser 
(or  hner)  hlter-width;  its  changes  are,  therefore,  less  dramatic  than  would  be  ex¬ 
pected  when  the  hlter-width  is  changed  by  a  factor  of  two.  No  such  memory  ehects 
are  present  in  the  Smagorinsky  model,  were  the  eddy  viscosity  is  more  tightly 
connected  to  the  hlter  width  Aj. 

This  is  clearly  observed  in  Figure  4.7,  which  shows  the  streamwise  development 
of  turbulent  kinetic  energy,  Ln,  and  the  normalized  eddy  viscosity.  We  observe 
sharp  discontinuities  in  the  CF-S  and  FC-S  cases  for  ut,  which  are  directly  related 
to  the  grid  size.  The  variable  hlter-width  is  in  this  case  not  as  benehcial,  somewhat 
smoothing  the  discontinuities  but  also  resulting  in  delayed  recovery  (especially  for 
the  CF-V  computation). 

The  diherence  with  the  single  grid  case  of  the  longitudinal  spectra  (Figure  4.8) 
is  more  signihcant  than  for  the  LDEV  case  as  well.  We  observe  a  signihcant  delay 
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Figure  4.6:  Streamwise  distribution  oi  (a)  turbulent  kinetic  energy  K]  (b)  integral 
length  scale  Ln;  and  (c)  normalized  eddy  viscosity,  vt/^-  LDEV  model  with  lower 
resolution,  48^  points  per  block  on  the  hne  side  and  24^  on  the  coarse  side. 
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Figure  4.7:  Streamwise  distribution  of  (a)  turbulent  kinetic  energy  K]  (b)  integral 
length  scale  Ln;  and  (c)  normalized  eddy  viscosity,  vt/v-  Smagorinsky  model. 


90 


Figure  4.8:  Longitudinal  spectra.  Smagorinsky  model.  Each  spectrum  is  shifted 
upward  by  100  units  for  clarity. 

in  the  regeneration  of  the  small-scale  energy  content  for  the  CF-V  case,  the  energy 
pileup  at  small  scales  for  the  FC-S  simulation  and  the  improved  behavior  of  the  flow 
immediately  after  the  discontinuity  in  the  FC-V  calculation. 

4.2.5  Two-level  calculations  with  explicit  filtering  of  the  advective 
term 

The  possible  beneht  of  hltering  explicitly  the  advective  term  in  the  momentum 
equations  derives  from  the  fact  that,  by  removing  a  band  of  the  highest  frequencies 
allowed  by  the  mesh,  the  truncation  and  aliasing  errors  can  be  reduced  [59].  For 
these  reasons  we  applied  explicit  hltering  together  with  the  LDEV  model  on  the 
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same  cases  discussed  before.  That  includes  two-level  calculations  where  the  grid 
transitions  from  96^  points  to  48^  points  in  the  Fine-to-Coarse  case  (FC)  or  vice- 
versa  in  the  Coarse-to-Fine  (CF)  situation.  In  order  to  obtain  similar  levels  of 
inflow  turbulent  kinetic  energy  as  in  Subsection  4.2.3,  the  dissipation  constant  for 
the  linearly  forced  isotropic  turbulence  runs  was  set  to  e  =  0.2U^/Lr.  The  resulting 
Reynolds  number  at  the  interface  location  was  Re\  ~  800.  The  same  variable- 
hlter  strategy  on  the  fine  side  discussed  before  was  employed  for  the  filtering  of  the 
advective  term. 

The  flow  visualizations  in  Figure  4.9  show  the  expected  trends  of  eddy  de¬ 
pletion  on  the  transition  to  coarser  mesh  for  the  sharp  FC,  and  scales  generation 
as  the  structures  are  advected  to  a  better  resolution  environment,  in  the  CF  cases. 
However,  we  also  notice  a  more  gradual  generation  of  small  scales  for  the  CF-S 
simulation.  Figure  4.9(a).  More  importantly,  for  the  FC-V  simulation  we  observe  a 
more  pronounced  increase  of  the  eddy  size  on  the  hne-grid  side,  compared  with  the 
case  with  no  explicit  hltering.  Figure  4.3(c).  The  scale  of  the  eddies  crossing  the 
interface  is,  thus,  increased;  the  better  resolution  of  these  eddies  results  in  better 
accuracy  as  the  flow  transitions  from  the  fine  to  the  coarse  grid. 

Consequently,  the  spectra  evolve  much  more  gradually  across  the  interface 
(Figure  4.10).  In  the  CF-S  case  the  high  wave-number  part  of  the  spectrum  is 
hlled  within  a  distance  of  0.7 (less  than  one  integral  scale  of  the  flow).  When 
the  hlter-width  is  decreased  gradually  (CF-V  simulation)  a  longer  transition  results, 
as  observed  in  all  previous  cases.  In  the  Fine-to-Coarse  calculations  we  observe 
no  energy  pileup  at  small  scale,  and  a  very  rapid  transition  for  the  discontinuous 
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Figure  4.9:  Isosurfaces  of  Q,  colored  by  the  value  of  the  eddy  viscosity.  LDEV  model 
with  explicit  hltering.  (a)  CF-S;  (b)  FC-S;  (c)  FC-V.  The  flow  is  from  top  right  to 
bottom  left. 
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Figure  4.10:  Longitudinal  spectra  for  LDEV  model  with  explicit  filtering.  Each 
spectrum  is  shifted  upward  by  1000  units  for  clarity. 

hlter-width  variation,  a  more  gradual  one  for  the  variable  hlter-width  case,  FC-V. 
Note  that  the  explicit  hltering  dampens  the  high-wavenumber  content  of  the  non¬ 
linear  term,  which  results  in  more  rapid  energy  decay  in  all  calculations  of  this  type 
(including  the  single-grid  ones). 

In  Figure  4.11  the  mean  turbulent  kinetic  energy,  integral  length,  and  eddy  vis¬ 
cosity  are  shown.  The  sharp  FC  transition,  as  in  previous  results,  exhibits  a  jump  as 
the  grid  transitions  from  hne  to  coarse;  interestingly,  however,  there  is  insignihcant 
energy  accumulation  on  the  fine  side  (Figure  4.10).  This  result  can  be  explained  by 
noting  that  explicit  hltering  produces  a  depletion  of  energy  on  the  structures  at  the 
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Figure  4.11:  Streamwise  distribution  of  (a)  turbulent  kinetic  energy  K]  (b)  integral 
length  scale  Ln;  (c)  normalized  eddy  viscosity,  LDEV  model  with  Explicit 

Filtering  of  advective  term. 
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grid  scale,  where  numerical  error  is  significant  (also  evident  in  Figure  4.10  at  the 
first  two  stations).  Since  the  smallest  scales  are  better  resolved  by  both  grid  levels, 
the  eddies  advected  across  the  interface  cause  less  of  a  perturbation  to  the  flow. 
The  longitudinal  spectra  in  Figure  4.10  show  that,  for  the  varying  filter  FC-V  case, 
a  gradual  reduction  of  energy  content  on  the  higher  wave-numbers  of  the  fine  grid 
is  achieved.  This  enables  a  natural  transition  to  the  natural  energy  distribution  for 
the  coarse  grid,  which  has  a  smaller  range  of  resolved  wave-numbers.  The  FC-S  case 
(consistent  with  the  previous  sections  results)  displays  a  discontinuous  transition  on 
the  spectra  as  the  resolved  wave-number  range  is  reduced.  The  coarse-to-fine  tran¬ 
sition  behaves  in  fashion  similar  to  that  seen  for  the  previous  LDEV  calculations, 
the  sharp  variation  being  superior  to  the  variable  filter-width  case,  yielding  a  faster 
recovery  of  fine  grid  statistical  data.  The  variable  filter  FC-V  calculation,  on  the 
other  hand  exhibits  a  faster  recovery  of  the  coarse  grid  statistics  than  the  sharp 
case,  for  both  integral  length-scale  and  turbulent  viscosity.  Figures  A.lllh)  and  (c). 

4.3  Flow  around  a  Sphere  at  Re= 10000 

In  this  section  we  apply  the  spatially  varying  LDEV  model  with  explicit  filter¬ 
ing  of  the  advective  term  to  LES  of  the  flow  around  a  sphere  at  Re  =  U^oD /v  =  10"^ 
{Uoo  is  the  freestream  velocity  and  D  is  the  diameter  of  the  sphere)  is  considered. 
This  Reynolds  number  falls  in  the  subcritical  regime,  where  laminar  boundary  layer 
separation  occurs  on  the  surface  of  the  sphere,  and  transition  happens  in  the  de¬ 
tached  shear  layers  in  the  wake.  We  have  seen  in  the  previous  Section  that  the 
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varying  filter  LDEV  model  with  explicit  filtering  of  convective  terms  is  suitable  for 
situations  where  turbulence  is  advected  across  derehnement  jumps.  That  is  the  case 
in  this  problem,  as  we  move  away  from  the  sphere  along  the  wake  region.  A  rectan¬ 
gular  domain  is  employed,  spanning  38D  in  the  streamwise  direction  and  24D  in  the 
cross  stream  directions.  The  sphere  is  centered  at  a  distance  of  lOD  downstream 
of  the  inflow  boundary.  The  inflow  velocity,  Uoo,  is  prescribed  at  this  boundary, 
and  a  convective  condition  used  at  the  outflow  boundary  [74].  Slip  wall  boundary 
conditions  are  used  in  the  other  two  directions. 

A  preliminary  computation  for  laminar  flow  around  a  sphere  at  Re  =  300  was 
done  on  a  similar  setting,  using  5  levels  of  rehnement  and  about  6.5  x  10®  points. 
The  resulting  grid  spacing  close  to  the  sphere  was  0.0195/1.  The  predicted  mean 
drag  coefficient  is  Cd  =  0.634,  and  its  oscillation  amplitude  C'^  =  0.0039.  These 
values  are  in  good  agreement  with  the  corresponding  results  from  computations  by 
Johnson  &  Patel  [47]  of  Cd  =  0.656  and  C'^  =  0.0035  respectively.  The  resulting 
Strouhal  number  is  St  =  0.132,  which  is  also  in  good  agreement  with  the  value  0.137 
reported  in  [47]. 

For  the  case  at  Re  =  10^  a  hner  grid  with  6  levels  of  rehnement  is  used  with 
11  X  7  X  7  blocks  at  level  0.  Each  block  consists  of  16^  cells.  The  grid  is  hnest 
near  the  surface  of  the  sphere  and  it  is  gradually  derehned  in  the  wake.  The  leaf- 
block  arrangement  at  an  x  —  z  plane  passing  through  the  center  of  the  sphere  is 
shown  in  Fig.  4.12.  The  cell  sizes  near  the  surface  of  the  sphere  are  approximately, 
Axi  =  0.003/1,  and  the  resulting  distance  of  the  hrst  grid  point  off  the  wall  is 
between  r~^  =  1.3  —  2.3  (following  [24],  we  dehne  r~^  =  rUrli'i  where  the  friction 
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Figure  4.12:  AMR  grid-blocks  for  flow  around  a  sphere  at  Re  =  10^.  An  x  —  z 
plane  at  y  =  0  is  shown.  Note  that  a  window  spanning  from  —5  <  a;  <  28  in  the 
x-direction  and  from  —5  <  <  5  in  the  z-direction  is  shown.  Six  levels  of  rehnement 

are  used  with  16^  points  per  block. 

velocity  is  assumed  to  be  Ur  =  0.04).  The  number  of  points  between  the  wall  and 
r~^  ~  10  is  between  5  —  8.  The  total  number  of  blocks  is  18400,  and  the  number  of 
leaf  blocks  16184  corresponding  to  66.3  million  points. 

The  equations  were  integrated  in  time  until  the  effect  of  initial  conditions  was 
eliminated.  Statistics  were  accumulated  over  two  shedding  cycles,  which  is  about 
lOUoo/D  eddy  turnover  times  or  20000  timesteps.  Although  the  sample  is  fairly 
small  it  is  sufficient  to  provide  reasonable  estimates  of  the  average  force  coefficients. 
LES  of  the  flow  around  a  sphere  at  the  same  Reynolds  number  performed  by  Con- 
stantinescu  et  al.  [24],  gave  a  mean  drag  coefficient  Cn  =  0.393.  The  same  value 
was  reported  in  the  LES  by  Yun  et  al.  [117].  Our  computation  gave  Cd  =  0.405, 
which  is  in  very  good  agreement  with  the  adobe  results.  In  Fig.  4.13  the  distribution 
of  average  pressure  coefficient,  Op,  and  skin  friction  coefficient,  O/,  are  shown  as 
a  function  of  the  azimuthal  angle,  0.  The  corresponding  results  from  the  LES  by 


Figure  4.13:  Variation  of  Cp  and  Cj  on  the  surface  of  the  sphere  as  a  function  of 
angle  0  from  front  stagnation  point.  Both  are  averaged  in  the  azimuthal  direction 

and  time.  - current  AMR-LES  calculation  at  Re  =  10^,  —  reference  LES  at 

Re  =  10*^  [24],  and  •  experiment  at  Re  =  1.6  x  10^  [1]. 

Const  ant  inescu  et  a/.  [24]  at  the  same  Reynolds  number,  and  the  experiments  by 
Achenbach  [1]  at  a  higher  subcritical  Reynolds  number.  Re  =  1.6  x  10^,  are  also 
shown  for  comparison.  In  general  the  agreement  is  good,  and  Cp  is  practically  the 
same  in  all  three  data  sets  for  0"  <  0  <  90°.  Both  simulations  also  agree  very 
well  for  90°  <  0  <  130°,  and  slightly  over-predict  Cp  reported  in  the  experiment. 
Further  behind  the  sphere  at,  130°  <  0  <  180°,  our  Cp  is  lower  by  is  approximately 
6%  of  the  peak  value,  probably  due  to  our  limited  statistical  sample. 

The  friction  coefficient,  Cf,  is  also  in  good  agreement  with  the  reference  data. 
In  this  case  discrepancies  among  the  different  datasets  are  larger,  due  to  larger 
errors  in  measuring  this  quantity  as  well  as  its  sensitivity  to  numerical  resolution. 
Our  mean  separation  occurs  at  0  =  91°  which  is  in  very  good  agreement  with  the 
value  of  0  =  90°  reported  in  [117].  Constantinescu  et  al.  [24]  is  reported  mean 
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Figure  4.14:  A  snapshot  of  the  instantaneous  flow  held  around  the  sphere  at  Re  = 
10'^.  (Top):  Vorticity,  oJy,  contours  at  an  x  —  2;  plane  through  the  center  of  the 
sphere.  Forty  contours  between  IDyD/Uoo  =  —20  (blue)  and  20  {red)  are  shown. 
(Bottom):  Q  isosurfaces  in  the  wake.  The  grid  distribution  is  also  shown  for  3  slices 
at  positions  x  =  0,  3D,  7D. 
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separation  at  0  =  85".  They  used  spherical  coordinates  and  a  grid  which  is  much 
hner  near  the  surface  of  the  sphere,  while  our  resolution  is  comparable  to  the  one 
utilized  by  Yun  et  al.  [117]. 

In  Fig.  4.14  an  instantaneous  snapshot  of  the  spanwise  vorticity,  uJy,  is  shown 
at  an  X  —  z  plane  passing  through  the  center  of  the  sphere.  The  shear  layer  insta¬ 
bility,  roll-up  and  subsequent  breakdown  can  be  observed.  The  smooth  variation 
of  the  vorticity  between  blocks  at  different  refinement  levels  should  also  be  noted, 
indicating  that  turbulent  eddies  are  unaffected  by  the  presence  of  block  boundaries. 
A  three-dimensional  view  of  such  eddies  is  also  shown  in  Fig.  4.14,  and  is  visual¬ 
ized  using  the  second  invariant  of  the  velocity  gradient  tensor,  Q.  The  AMR  block 
structure  at  three  different  planes  in  the  wake  region  [x  =  0,  x  =  3D  and  x  =  7D) 
is  also  shown.  As  the  rehnement  level  is  decreased  for  increasing  x,  the  reduction 
in  the  range  of  scales  resolved  by  the  grid  is  evident. 
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Chapter  5 


Applications  to  flapping  flight 

Historically,  flapping  flight  has  received  great  attention  from  the  biology,  engi¬ 
neering  and  science  disciplines.  For  biologists,  the  interest  relies  on  the  mechanisms 
by  which  biological  motor,  sensor  and  control  has  evolved  and  is  employed  by  flying 
insects  in  everyday  survival  [87].  Engineering  research  on  flying  has  been  focused 
hxed  wing  aerodynamics,  where  their  detailed  understanding  and  control  has  driven 
research  over  the  past  century,  mainly  as  a  quest  to  serve  the  developing  aeronautical 
and  aerospace  industries.  However,  advances  in  material  technologies  and  actuator 
miniaturization  in  the  past  decade  has  sparked  new  interest  in  small  scale  flapping 
wing  devices  [91].  As  the  size  of  the  vehicle  and  characteristic  Reynolds  number 
{Re)  of  the  flow  decrease,  friction  forces  are  enhanced,  and  fixed  wing  aerodynamics 
fails  to  account  for  the  forces  required  to  sustain  flapping  flight  and  maneuvering  in 
an  efficient  manner  [7]. 

Several  experimental  and  computational  investigations  have  been  carried  out 
over  the  past  years,  to  understand  the  unsteady  mechanisms  that  allow  the  genera¬ 
tion  of  enough  aerodynamic  force  to  hover  and  maneuver.  Experimental  studies  have 
given  insights  on  the  kinematic  and  kinetic  aspects  of  flapping  flight  on  thetered  and 
free  flying  animals  ([108],  [38]),  as  well  as  dynamically  scaled  robotic  devices  ([35], 
[27],  [88]).  Viscous  incompressible  flow  simulations  of  flapping  wing  models  has  also 
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been  reported  ([58],  [93],  [82],  [105],  [106],  [107],  [83]).  These  studies  have  been  pri¬ 
marily  used  as  an  extension  of  experiments  to  obtain  a  more  detailed  view  of  the  flow 
helds  and  the  large  vortices  that  dominate  the  production  of  aerodynamics  forces. 
Liu  and  Kawachi  [58],  for  example,  conducted  one  of  the  hrst  computations  of  flow 
around  a  flapping  Manduca  Sexta  moth  wing  at  Re  =  4000.  Due  to  the  lack  of  nu¬ 
merical  resolution  only  some  of  the  large  scale  structures  had  some  similarity  to  the 
experimental  data  in  [109].  Ramamurti  and  Sandberg  performed  similar  study  for 
the  smaller  Drosophila  melanogaster  fly  at  Re  =  136.  Most  other  computations  in 
the  literature  utilize  two-dimensional  (2D)  models  ([105],  [106]).  In  addition,  insect 
wings  are  complex  flexible  structures  that  add  signihcant  complexities  to  numerical 
modeling. 

Today  most  of  the  computational  studies  on  flapping  flight  focus  on  the  fluid 
dynamics  of  rigid  flapping  wings  with  prescribed  kinematics,  or  examine  the  wing  as 
a  structure,  modeling  the  fluid  loads  with  added-mass  and  damping  parameters.  The 
coupling  between  the  fluid  and  the  structure  is  either  neglected  or  taken  into  account 
by  extremely  coarse  models.  In  general,  a  flexible  wing  structure  deforms  under  the 
influence  of  aerodynamic  forces,  elastic  deformation  forces,  and  also  inertia  forces 
due  to  the  accelerations  of  its  own  mass  ([33],  [36],  [53],  [23]).  The  deformation 
is  given,  to  a  large  extent,  by  the  internal  distribution  of  compliant  components 
and  mechanisms  ([110]).  The  relative  contributions  of  aerodynamic  and  inertial- 
elastic  forces  on  wing  deformation  of  the  Manduca  sexta  hawkmoth,  for  example, 
were  assessed  by  Combes  and  Daniel  (2003)  [23].  They  concluded  that,  for  Manduca 
sexta  the  equilibrium  of  forces  during  the  wing’s  motion  is  mainly  due  to  inertia  and 
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elastic  forces,  and  the  aerodynamic  loads  had  a  minor  role  in  the  form  of  damping. 
For  this  insect  the  typical  ratio  of  wing  inertial  to  aerodynamic  forces  in  hovering 
flight  is  around  7.  In  other  species  this  parameter  is  closer  to  one  [36],  and  more 
complex  aeroelastic  interactions  occur.  Strongly  coupled  interactions  between  fluid 
and  wing  dynamics  require,  accurate  and  robust  accurate  FSI  solvers  [114]  such  as 
the  one  outlined  in  Chapter  3. 

In  this  chapter  we  will  utilize  the  AMR/FSI  solver  to  study  two  important 
aspects  of  flapping  flight:  i)  the  effects  of  flexibility  on  the  dynamics  of  hovering 
flexible  prohle,  and  in  particular  the  nonlinearities  of  the  structure  and  viscous 
flow  play  a  part  in  the  evolution  of  the  system,  ii)  the  dynamics  of  free  longitudinal 
flight  in  a  realistic  three-dimensional  conhguration.  Details  are  given  in  the  following 
sections. 

5.1  Influence  of  Flexibility  on  the  performance  of  two-dimensional 
flapping  flexible  airfoils 

In  this  section  we  will  study  the  the  effects  of  structural  flexibility  on  the 
aerodynamic  performance  of  a  flapping  wing  for  a  range  of  Reynolds  numbers.  We 
will  perform  two-dimensional  computations  in  order  to  explore  a  wide  range  of 
parameters  in  a  cost-efficient  manner.  In  particular,  a  representative  section  of  a 
three-dimensional  wing  (two-dimensional  foil)  is  considered,  and  spanwise  bending 
and  torsion  flexibility  are  neglected.  Such  a  structure  can  be  represented  by  two 
plates  connected  with  a  torsion  spring  to  account  for  deformation  in  the  chordwise 
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direction  (see  Figure  5.1a).  This  system  has  four  degrees  of  freedom,  which  are 
effectively  reduced  to  one  by  prescribing  harmonic  hovering  motions  of  one  of  the 
plates.  The  plates  are  rigid  and  the  large  angular  deformations  give  rise  to  cubic  and 
higher  order  odd  non-linearities  in  the  governing  equations,  similar  to  those  seen  in 
a  pendulum  and  in  flexible  beams  [5].  In  a  sense,  one  could  consider  the  two  plates 
as  a  double  pendulum  with  a  torsion  spring. 

5.1.1  The  structural  model 

A  schematic  of  the  structural  system  is  shown  in  Figure  5.1a.  It  consists  of 
two  rigid  links  A  and  B,  which  are  joined  at  the  center,  b,  by  a  pin  that  contains 
a  linear  torsion  spring.  In  this  model,  flexibility  is  concentrated  at  one  discrete 
location  of  the  system,  and  inclusion  of  elastic  links  will  allow  chordwise  variations 
of  stiffness  and  mass  to  be  accounted  for.  The  two  links  are  covered  by  a  set  of 
aerodynamic  surfaces  that  dehne  the  boundary  between  the  airfoil  and  the  fluid,  and 
deform  as  the  angle  between  them  changes.  The  aerodynamic  surfaces  consist  of  two 
rigid  segments,  Rsa  and  Rsb  (see  Figure  5.1b),  and  two  segments  that  dynamically 
deform  according  to  the  angle  between  the  two  links.  The  deformation  is  prescribed 
by  htting  the  Hermite  interpolation  polynomials  cl-c2  and  c3-c4.  We  have  found 
that  this  modeling  is  robust  and  helps  maintain  the  smoothness  of  the  surface  even 
for  large  values  of  the  angle  between  the  plates. 

The  dimensionless  form  of  the  equations  governing  the  motion  of  the  structural 
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(b) 


Figure  5.1:  (a)  Two-link  profile.  The  rigid  links  A  and  B  are  connected  at  hinge 
&  by  a  torsional  spring  of  stiffness  k.  The  variables  x(t),  y(t),  6{t)  and  a{t)  are 
the  generalized  coordinates  used  to  describe  the  wings  motion.  In  the  hovering 
simulations,  x(t),  y{t)  and  6{t)  are  prescribed  and  a{t)  is  left  as  only  degree  of 
freedom  required  to  dehne  the  system,  (b)  Decomposition  of  the  prohles  surfaces 
into  rigid  and  deformable  sections.  The  two  rigid  surfaces  Rsa  and  Rsh  are  connected 
at  points  Ci,  C2,  C3  and  C4  by  Hermite  interpolating  polynomials  Hgi  and  Hs2-  Also, 
rriA,  rns  are  the  total  masses  of  links  A  and  S,  and  Vb  are  the  respective 
distances  from  their  centers  of  mass  to  junction  b. 
system  shown  in  Fig.  5.1a  can  be  written  as: 
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(5.1) 


where  x{t),  y(t),  6{t)  are,  respectively,  the  joint  horizontal  displacement,  joint  verti¬ 
cal  displacement,  and  orientation  angle  of  link  B  measured  from  the  inertial  reference 
frame  N  as  shown  in  Figure  5.1a,  and  a{t)  is  the  deflection  angle  between  links  A 
and  B.  Here,  rrij  is  the  total  mass  of  plate  i  {i  =  A,  B)]  r]i  is  the  distance  from  the 
pin  joint  to  the  center  of  mass  of  bar  i]  is  the  moment  of  inertia  of  link  i  respect 
to  hinge  point  b.  Also,  Qx,  Qy  are  the  fluid  forces  in  the  x  and  y  directions  re¬ 
spectively,  and  Qe,  Qa  are  the  fluid  moments  acting  on  the  generalized  coordinates 
9(t)  and  a(t).  The  quantities  gx,  Qy,  ge,  ga  are  the  corresponding  contributions  of 
centrifugal,  elastic  and  gravity  forces.  The  equations  are  made  dimensionless  by 
using  the  chord  length  of  the  undeformed  airfoil,  as  the  reference  length  scale, 
and  the  maximum  translational  velocity  of  the  junction  b  between  the  two  links,  Uc 
as  the  reference  velocity.  The  Reynolds  number  is  dehned  as  Re  =  LJJd v-,  where  v 
is  the  fluid  kinematic  viscosity.  The  fluid  forces  and  moments  are  determined  from 
equations  (2.1). 

In  the  present  numerical  experiments,  the  translational  motions  of  the  junction 
as  well  as  the  orientation  of  link  B  are  prescribed.  With  these  prescribed  motions, 
the  four  degrees  of  freedom  of  the  system  can  be  effectively  reduced  to  the  deflection 
angle  a{t)  between  plates  A  and  B.  Thus,  the  overall  deformation  of  the  wing  section 
is  determined  by  the  deflection  angle  a{t),  which  is  governed  by  the  following  reduced 
form  of  equations  (5.1): 

Iach  +  ka  =  —Ia9  +  mAgAsin  {9  +  a)  x  +  Qa  (5.2) 

Equation  (5.2)  resembles  the  equation  governing  a  harmonic  oscillator  with  forcing 
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due  to  the  prescribed  kinematics  and  the  fluid  forces  (e.g.  reference  [71]).  The 
non-linearities  arise  from  the  sin{d  +  a)  term  due  to  the  kinematics  and  the  fluid 
forcing.  For  this  study,  we  only  take  into  account  the  fluid  damping  which  arises 
through  the  fluid  moment  Qa-  It  should  be  noted  that  selecting  a  proper  structural 
damping  model  is  far  from  trivial,  and  this  is  an  active  research  topic  in  structural 
biomechanics.  Damping  models  for  insect  wings  are  relatively  few  [e.g.  classical 
viscous  damping  model  used  by  Herbert  [22]  and  the  viscoelastic  model  used  by 
Bao  and  colleagues  [13]]  and  the  existing  models  require  a  fair  amount  of  empirical 
information. 

5.1.2  Parametric  space  and  computational  setup 

To  prescribe  the  translational  motions  of  the  junction  and  the  orientation  of 
link  B,  we  define  the  states  x(t),  y(t)  and  9(t)  as: 

6{t)  =  6*0  -|-  ^1  —  j  'ysin{uft  +  0),  (5.3) 

where  Aq  is  the  stroke  length  of  the  pin  point,  6q  is  the  mean  orientation  angle  for 
link  H,  7  is  the  rotation  amplitude,  cu/  is  the  frequency  of  the  prescribing  or  forcing 
oscillation  and  0  is  the  phase  angle  between  x{t)  and  9{t).  The  exponential  terms 
were  used  in  order  to  reduce  transient  effects  [23] .  The  time  constant  was  chosen  as 
T  =  I.Gtt/o;/  because  99.8%  of  the  prescribed  amplitude  was  reached  after  a  time 
length  of  5  periods.  The  following  parameters  corresponding  to  symmetric  hovering 
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were  selected  [107]: 


7lo  =  2.8;  ^o  =  -|;  7  =  |;  0  =  0  (5.4) 

Based  on  the  adopted  normalization,  the  problem  is  completely  dehned  by 
the  density  ratio  Pb/Pf,  the  frequency  ratio  Uf/un  and  the  Reynolds  number  Re. 
Here,  pb  is  the  density  of  the  wings  material  and  ojn  =  a/ {k/ Ia)  is  the  linear  natural 
frequency  of  the  oscillator  (equation  (5.2)).  The  frequency  ratio  ojf/ojn  is  used  to 
characterize  the  flexibility  of  the  wing  section. 

Three  Reynolds  numbers  were  considered  {Re  =  75,  250  and  1000)  to  inves¬ 
tigate  the  effect  of  the  reduction  in  viscous  dissipation  on  the  system  dynamics. 
The  mass  ratio  was  set  to  Pa/pj  =  25,  as  this  value  provided  a  ratio  close  to  2 
for  the  maximum  translational  inertial  force  over  maximum  drag  force  at  Re  =  75 
for  the  chosen  geometry  and  kinematics.  The  above  ratio  was  determined  through 
numerical  experiments  with  the  rigid  wing.  To  compute  the  maximum  horizontal 
translational  inertial  force,  the  total  wing  mass  was  multiplied  by  the  maximum 
acceleration  determined  from  the  second  derivative  of  x{t)  in  equation  (5.3).  The 
value  of  peak  drag  force,  on  the  other  hand,  was  obtained  from  the  rigid  wing  simu¬ 
lation  at  Re  =  75.  The  wing  has  a  thickness  of  10%  of  the  undeformed  chord  length 
and  circularly  formed  edges. 

The  frequency  ratio  in  the  case  of  Re  =  75  was  set  to  uif/uin  =  1/2,  1/2.5, 
1/3,  1/3.5,  1/4,  1/6,  and  for  Re  =  250,  1000  was  set  to  1/2,  1/3,  1/4,  1/6.  The 
resulting  range  of  maximum  deflection  angles  varied  from  10°  up  to  70°.  Also,  the 
rigid  prohle  problem  (no  angular  deformation  between  the  links)  was  run  for  each  of 
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the  Reynolds  numbers  mentioned.  We  note  that  for  frequency  ratios  below  1/2,  the 
computations  would  fail  since  the  two  plates  collide  during  rotation.  This  limitation 
arises  from  the  fact  that  the  flexibility  in  the  present  model  is  concentrated  at  the 
hinge  point  and  the  distributed  chordwise  variations  of  stiffness  and  mass  are  not 
accounted  for.  A  flexible  beam  model  and/or  inclusion  of  structural  damping  may 
help  to  address  this  issue  and  enable  computations  with  frequency  ratios  of  about 
one. 

The  grids  adopted  for  the  Re  =  75  simulations  were  tailored  to  maintain  a 
high  rehnement  zone  in  the  areas  where  boundary  and  detached  shear  layers,  as  well 
as  the  wake  vortical  structures  are  present.  The  rigid  wing  was  set  to  move  in  the 
center  of  a  box  with  dimensions  30Lc  x  30Lc,  in  order  to  minimize  interference  effects 
from  the  domain  boundaries.  For  all  this  simulations  we  chose  to  use  a  single  block 
stretched  grid,  where  a  fast  direct  Poisson  solver  could  be  used  [95] .  A  near-uniform 
grid  zone  was  generated  near  the  center,  where  the  motions  of  the  two-link  system 
took  place,  and  this  zone  was  stretched  towards  the  boundaries.  The  stretching 
towards  the  boundaries  of  the  domain  was  obtained  using  the  hyperbolic  tangent 
function  in  each  direction. 

For  the  Re  =  75  simulations,  the  uniform  grid  zone  had  a  local  cell  size 
of  Ax  =  Ay  =  0.0038Lc,  and  the  total  number  of  points  was  1229  x  551  along 
the  X  and  y  directions,  respectively.  Through  grid  rehnement  studies,  we  found 
that  the  above  resolution  was  sufficient  to  capture  all  how  features.  In  Figure  5.2, 
computationally  obtained  aerodynamic  forces  are  shown  for  approximately  half  the 
resolution  throughout  the  computational  domain  (total  number  of  points  was  664  x 
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400)  and  the  same  forcing  conditions  (i.e.  r  =  0  in  equations  (5.3))  and  Reynolds 
number  as  that  for  the  baseline  rigid  wing  computation.  The  corresponding  lift 
and  drag  coefficients  determined  in  the  computations  of  Wang  and  colleagues  [107], 
where  a  hovering  ellipse  with  the  same  kinematics  is  considered  instead  of  a  plate, 
are  also  included  in  the  hgure.  The  agreement  between  the  results  obtained  with  the 
two  different  grids  is  good,  with  a  maximum  difference  of  around  3%.  Despite  the 
differences  in  the  wing-section  shapes,  after  the  initial  transients  (t/T  >  2),  good 
agreement  is  also  seen  with  the  results  obtained  by  Wang  and  colleagues  [107]. 

For  the  results  of  Figure  5.2,  within  the  boundary  layers  on  the  link  surface, 
we  estimated  the  number  of  grid  points  to  be  approximately  8  and  16  for  the  coarse 
and  hne  grids,  respectively.  As  the  Reynolds  number  increases,  the  boundary  layer 
thickness  is  expected  to  decrease  as  a/I/Rc,  and  in  order  to  keep  the  resolution 
within  the  above  range,  a  grid  with  1320  x  1038  points  was  found  to  be  sufficient 
for  both  Re  =  250  and  Re  =  1000  cases.  All  the  results  presented  in  this  article 
were  obtained  with  a  1229  x  551  grid  for  the  Re  =  75  simulations  and  a  1320  x  1038 
grid  for  the  Re  =  250  and  Re  =  1000  simulations.  The  governing  equations  were 
integrated  for  a  time  length  of  14  periods,  21  periods  and  15  periods,  for  Re  =  75, 
250  and  1000,  respectively.  The  time-averaged  quantities  were  computed  over  the 
last  7  periods,  13  periods  and  10  periods,  for  Re  =  75,  250  and  1000,  respectively. 

5.1.3  Results 

Aerodynamic  quantities 
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Figure  5.2:  Time  histories  of  lift  and  drag  force  coefficients  (C^,  Cd)  for  a  symmetric 
harmonic  hovering  rigid  link  at  Re  =  75  and  two  different  grid  resolutions.  Blue 
line,  rigid  link,  embedded  boundary  grid  1229  x  551;  green  line,  embedded  boundary 
grid  666  x  402;  and  red  line,  data  from  Wang  and  colleagues  [107].  t,  time;  T  is  the 
prescribed  motion  time  period. 

Given  the  complexity  of  the  flexible  prohle  conhguration  as  a  function  of  defor¬ 
mation,  a  simple,  yet  consistent,  normalization  strategy  is  adopted  to  obtain  force 
and  power  nondimensional  coefficients.  That  is,  the  lift  and  drag  coefficients  are 
dehned  by: 
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where  Ql(t)  and  Q^it)  are  dimensional  quantities  and  Qx(t)  and  Qy(t)  are  non- 
dimensional  quantities.  Once  the  equation  for  the  deformation  a{t)  is  solved  at 
each  integration  step,  the  driving  forces  in  the  prescribed  generalized  coordinates 
x(t),  y(t)  and  9(t)  are  computed  from  5.1. 

The  total  power  input,  sum  of  horizontal  translational  power  and  rotational 
power,  is  computed  from 

Ptr  =  Rx{t)  X  x{t) 

Prot  =  Re(t)  X  9(t),  (5.6) 

where  Ptr  and  Prot  are  the  translational  and  rotational  power  inputs  at  the  hinge 
b,  Rx{t)  is  the  driving  force  in  the  x  direction,  and  Relt)  is  the  driving  moment  in 
the  6*(t)  angular  direction.  In  an  ideal  case  where  the  driving  mechanism  on  the  pin 
is  perfectly  elastic,  all  the  negative  power  exerted  on  the  mechanism  can  be  stored 
as  potential  deformation  energy  for  later  use,  adding  to  the  profiles  aerodynamic 
efficiency.  Here,  following  Berman  and  Wang  [21],  a  conservative  approach  is  used. 
We  suppose  that  negative  power  is  not  reusable  and  will  not  be  taken  into  account 
in  the  computation  of  time  averages.  Then; 

I  Ptr,  if  Ptr>0  I  Prot,  if  Prot  >  0 

Ptr=<  ,  Prot  =  <  (5.7) 

0,  if  Ptr  <0  ^  0,  if  Prot  <  0 

The  power  coefficient  is  dehned  as: 

Cpw(t)  =  (^<y)  +  =  2(P,Jt)  +  P„(t))  (5,8) 

2-  PfUf-  Lo 

In  Figure  5.3  the  dependence  of  mean  values  of  Cl,  Co,  the  ratio  Cl/ Co  and 


113 


Cl/Cpw  with  respect  to  the  frequency  parameter  Uf/un  are  shown.  The  rigid  wing 
has  zero  torsion  stiffness  equivalent  to  =  0. 

For  all  Reynolds  numbers  considered,  lift  and  drag  coefficients  exhibit  a  peak 
at  the  value  of  ujfjujn  =  1/3.  The  performance  parameter  CL/CD  also  exhibits  a 
prominent  peak  at  this  frequency  ratio.  For  Re  =  75,  250  and  1000  increases  of 
28%,  23%  and  21%  respect  to  the  rigid  prohle  performance  were  observed. 

The  variations  of  the  aerodynamic  quantities  with  respect  to  the  frequency 
ratio  show  similar  characteristics  for  all  three  Reynolds  numbers. The  decrease  in 
viscous  forces  produces  a  decrease  in  drag  and  an  increase  in  lift  over  the  frequency 
parameter  range,  therefore,  Ci/Cp  is  generally  higher  for  the  higher  Re  cases.  It 
is  interesting  to  note  the  striking  difference  between  the  graph  of  Cl/Cd  obtained 
for  the  Re  =  75  case  and  those  obtained  for  the  higher  Reynolds  numbers.  For  the 
lowest  Reynolds  nnmber  and  ujflujn  =  1/4,  the  above  ratio  is  over  13%  higher  than 
that  obtained  for  the  rigid  wing,  while  for  Re  =  250  it  is  only  increased  by  0.5%. 
In  Figure  5.3c,  it  can  be  seen  that  for  Uf/uin  =  1/3,  the  performance  ratio  Ci/Cpw 
is  39%  and  28%  higher  than  that  obtained  for  the  rigid  wing  for  theRe  =  75  and 
Re  =  250  cases,  respectively.  Interestingly,  this  measure  is  only  abont  13%  higher 
than  that  obtained  for  the  rigid  case  at  Re  =  1000. 
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Figure  5.3:  (a)  Mean  values  of  Cl  and  Cd  as  a  function  of  the  frequency  ratio 
ujf/ujn-  blue  circle,  Cl  at  Re  =  75;  red  circle,  Cl  at  Re  =  250;  black  circle.  Cl 
at  Re  =  1000;  blue  diamond,  Cd  at  Re  =  75;  red  diamond,  Cd  at  Re  =  250; 
and  black  diamond,  Cd  at  Re  =  1000.  (h)  Ratio  of  mean  Cl/Cd  versus  ujf/uin'- 
blue  circle.  Re  =  75;  red  diamond.  Re  =  250;  and  black  triangle.  Re  =  1000.  (c) 
Mean  lift  coefficient  per  unit  of  driving  power  coefficient  {Cpw)  versus  ujf/ujn'-  same 
definitions  as  in  (b). 

In  Figure  5.4,  the  time  histories  of  the  lift  and  drag  coefficients  are  shown  for 
all  Reynolds  numbers  at  three  selected  frequency  ratios.  The  effects  of  ffexibility 
are  noticeable  on  the  lift-force  peaks  at  the  initiation  of  the  stroke  (indicated  with 
a  black  arrow  in  Fig.  5.4a).  For  Re  =  75  and  ujf/uJn  =  1/2,  corresponding  to  the 
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Figure  5.4:  Time  histories  of  lift  and  drag  coefficients  for  Re  =  75,  Re  =  250  and 
Re  =  1000.  (a)  lift  coefficient  and  (b)  drag  coefficient;  blue  dashed  line,  rigid  wing; 
red  dashed  line,  flexible  wing  with  ujf/ujn  =  1/2;  green  line,  flexible  wing  with 
uif  /uJn  =  1/3;  black  dash-dot  line,  flexible  wing  with  ujf  /ujn  =  1/4. 

most  flexible  airfoil,  this  peak  is  negative,  while  for  the  rigid  case  the  coefficient 
of  lift  peaks  at  0.5.  For  all  cases  in  between,  the  enhancement  of  the  mean  lift 
force  seen  in  Figure  5.3a  comes  from  the  gradual  increase  of  this  peak,  which  is 
at  0.83  and  1.28  for  uif/ujn  =  1/4  and  oJf/oJn  =  1/3,  respectively.  For  the  latter 
frequency  ratio,  where  a  structural  non-linear  resonance  occurs,  the  lift  peak  is  also 
delayed  and  nearly  coincides  with  the  translational  lift  peak.  This  is  translated  into 
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a  larger  area  under  the  lift  curve  per  period  and  a  larger  time-averaged  Cl  value. 
The  temporal  variations  of  the  lift  and  drag  coefficients  for  Re  =  250  and  1000  are 
more  complex  than  the  variations  seen  at  Re  =  75,  and  the  periodicity  is  practically 
lost.  Still,  in  an  average  sense,  negative  lift  peaks  after  stroke  reversal  and  larger 
translational  lift  peaks  are  seen  when  Ufjun  =  1/2.  Also,  a  widened  two-peak  lift 
cnrve  is  observed  when  Uf/un  =  1/3. 

Flow  strnctnres 

In  order  to  relate  the  temporal  variations  of  the  lift  and  drag  forces  to  spe- 
cihc  flow  strnctnres,  we  examine  several  realizations  of  the  instantaneous  flow  helds. 
In  Figure  5.5,  vorticity  isolines  are  shown  for  the  rigid  and  Uf/un  =  1/3  cases  at 
Re  =  250.  For  clarity,  the  lift  coefficient  variation  has  been  added  (Figure  5.5k), 
together  with  the  temporal  variation  of  the  phase-averaged  circulation  of  the  most 
important  vortical  strnctnres  generated  dnring  a  flapping  cycle.  These  are  the  lead¬ 
ing  edge  vortex  (LEV)  shown  in  Figure  5.5a,  the  end  of  stroke  vortex  (ESV)  shown 
in  Figure  5.5c  and  trailing  edge  vortex  (TEV)  shown  in  Figure  5.5e.  The  circulation 
of  each  of  these  vortices  was  compnted  as  a  fnnction  of  time  by  direct  integration  of 
the  vorticity  within  a  given  threshold  contonr  aronnd  each  vortex.  The  selection  of 
the  threshold  contonr,  althongh  arbitrary,  was  consistently  taken  to  be  the  lowest 
closed  vorticity  isoline  in  the  vicinity  of  the  given  vortex.  The  range  of  contours  was 
given  by  48  isolines  within  an  interval  of  10  vorticity  units  placed  around  a  hand 
picked  value  at  each  time.  This  value  was  obtained  on  a  position  of  the  shear  layer 
associated  with  the  given  vortex  and  close  to  it.  Then,  the  circnlation  of  each  vortex 
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was  first  computed  in  time  for  both  forward  and  reverse  strokes  over  nine  periods, 
and  then,  both  of  these  quantities  were  averaged.  The  number  of  frames  available 
per  period  for  the  simulations  was  40. 

As  the  flexible  wing  approaches  the  end  of  the  stroke  in  Figure  5.5a,  it  exhibits 
different  rotation  velocities  on  the  two  components  A  and  B.  The  driven  link  B 
rotates  with  the  prescribed  angular  velocity  0{t),  and  the  lower  link  A  rotates  with 
an  angular  speed  [9(t)+  ait)].  The  added  angular  speed  d(t)  affects  the  overall 
dynamics  at  stroke  reversal.  First,  the  camber  generated  by  the  angular  deformation 
a{t)  at  the  end  of  the  stroke  (see  Figure  5.5b)  reorients  the  zero  lift  direction  on  the 
wing  and  enhances  wake  capture  effects.  This  enhancement  mechanism  is  analogous 
to  the  one  produced  by  orientation  advancement  in  rigid  wings  [107]  .  It  is  important 
to  note  that  an  excessive  degree  of  flexibility  (low  frequency  ratios)  produces  a  large 
camber  at  stroke  reversal,  which  has  a  negative  effect  on  the  lift  production  (e.g.  at 
Ufjun  =  1/3  in  Figure  5.4).  The  evolution  and  strength  of  the  LEV  on  the  other 
hand  (see  Figure  5.5a)  is  only  a  weak  function  of  the  wings  flexibility.  The  formation 
time  as  well  as  the  maximum  circulation  shown  in  Figure  5.51  are  approximately  the 
same  for  the  rigid  and  the  flexible  wings. 

Another  effect  of  the  higher  rotation  speeds  at  the  trailing  edge  for  the  flexible 
wings  is  the  formation  of  a  stronger  shear  layer,  where  higher  vorticity  values  are 
found.  This  shear  layer  rolls  up  into  a  stronger  ESV  as  can  be  seen  in  Figure  5.5c. 
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Figure  5.5:  Behavior  of  rigid  and  flexible  profile  with  =  1/3  at  Re  =  75  during 

stroke  reversal,  (a-e)  Vorticity  contours  for  flexible  prohle  at  t/T  =  —0.0491,  0.0009, 
0.0759,  0.1760,  and  0.2260;  80  contours,  cUmin  =  —10,  (Umav  =  10.  (f-j)  Vorticity 
contours  for  rigid  prohle  at  same  time  locations.  White  dashed  lines  denote  end 
of  stroke  position,  (k)  lift  coefficient  history.  (1)  Circulation  histories  of  leading 
edge  vortex  (LEV),  end  of  stroke  vortex  (ESV),  and  trailing  edge  vortex  (TEV).  F, 
circulation. 
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On  examining  the  ESV  circulation  plots  (Figure  5.5k),  one  finds  that  the 
strength  and  life  span  are  significantly  enhanced  when  compared  with  those  of  a 
rigid  wing.  The  ESV  pinches  off  later,  forming  a  pair  of  counter-rotating  vortices 
together  with  the  LEV.  This  vortex  pair  generates  flow  directed  towards  the  wing, 
enhancing  the  wake-capturing  effects.  This  is  more  clearly  reflected  in  the  lift  coef¬ 
ficient  evolution  shown  in  Figure  5.5k. 

In  contrast  to  the  rigid  wing,  where  the  lift  curve  reaches  a  maximum  (point  H) 
and  starts  to  decrease,  for  the  flexible  wing  the  production  of  increased  lift  continues 
for  longer  (point  C).  Once  a  flexible  wings  deflection  has  reached  its  maximum, 
the  elastic  energy  stored  in  the  torsion  spring  is  released  to  generate  a  restoring 
motion,  whose  timing  again  depends  on  the  degree  of  flexibility  of  the  structure 
(see  Figures  5.5d  and  5.5e).  This  restoring  motion  produces  a  dynamical  change 
in  the  wings  camber  with  a  resulting  increase  in  the  fluid  forces.  This  also  affects 
the  formation  and  growth  of  the  TEV.  It  is  well  established  that  this  flow  structure 
generates  a  low-pressure  zone,  which  translates  into  increasing  forces  up  to  the 
pinch-off  time  (see,  for  example,  reference  [106]).  It  can  be  seen  in  Figures  5.5k-l, 
that  the  time  at  which  the  TEV  pinches  off  is  correlated  with  the  translational  force 
peak;  this  peak  is  advanced  in  the  uif/un  =  1/3  case  when  compared  with  that  for 
the  rigid  wing. 

In  Figure  5.6,  a  quantitative  comparison  of  the  LEV,  ESV  and  TEV  dynamics 
is  shown  for  different  flexibilities  at  Re=75  by  determining  their  average  circulations 
as  a  function  of  time. 
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Figure  5.6:  Averaged  circulations  as  a  function  of  time  with  respect  to  stroke  reversal 
at  Re  =  75;  blue  line,  TEV;  green  dashed  line,  LEV;  and  red  dashed-dot  line,  ESV. 
(a)  Rigid  airfoil,  (b)  Flexible  wing  with  ujf/ujn  =  1/2,  (c)  Flexible  airfoil  with 
ojf/ojn  =  1/3,  and  (d)  Flexible  profile  with  ojf/ujn  =  1/4. 

The  maximum  averaged  circulation  for  each  vortex  and  the  time  at  which  it 
occurs  with  respect  to  the  stroke  reversal  are  provided  in  Table  5.1.  As  expected, 
from  what  was  observed  in  Figure  5.5,  the  LEV  dynamics  is  similar  for  all  frequency 
ratios  in  terms  of  both  strength  and  timing.  The  TEV  on  the  other  hand,  attains  a 
higher  maximum  circulation  as  the  wing  becomes  more  flexible.  However,  the  time  it 
takes  to  reach  this  maximum  circulation  is  shortest  for  ujf/ujn  =  1/3,  where  the  best 
aerodynamic  performance  is  seen.  For  the  ESV  vortex,  the  maximum  circulation 
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^  f  / 

r  TEV 

Time  TEV 

r  LEV 

Time  LEV 

r  ESV 

Time  ESV 

Rigid 

2.44 

0.78 

1.30 

-0.02 

0.44 

0.001 

1/2 

3.14 

0.75 

1.30 

0.03 

0.31 

0.18 

1/3 

2.78 

0.58 

1.25 

0.001 

0.72 

0.10 

1/4 

2.29 

0.68 

1.36 

0.001 

0.87 

0.03 

1/6 

2.36 

0.80 

1.43 

0.001 

0.59 

0.02 

Table  5.1:  Values  of  maximum  mean  circulation  T  obtained  for  TEV,  LEV,  and  ESV 
at  Re  =  75.  The  circulations  were  obtained  by  integrating  the  vorticity  function  for 
each  vortex  for  two  strokes  and  taking  the  average  value  for  each  time  frame.  Time 
is  dehned  as  in  Figure  5.6  and  corresponds  to  the  maximum  mean  circulation. 

increases  for  the  frequency  ratios  ujfjujn  =  1/3  and  1/4.  The  peak  circulation  for 
f/n=l/3  is  20%  lower  than  that  obtained  for  f/n=l/4,  but  the  deflection  at  stroke 
reversal  in  terms  of  the  maximum  deformation  angle  is  90%  larger  when  cnj/cn^  =  1/3 
(56  deg.  ioi  ojf/ojn  =  1/3  and  29  deg.  for  Uf/un  =  1/4).  This  is  translated  into  a 
larger  projected  area  contribnting  to  the  lift  force.  Also,  as  seen  from  Figure  5.6, 
the  time  delay  of  the  peak  circulation  of  the  ESV  vortex  is  increased  as  the  wing 
becomes  more  flexible. 

A  more  direct  illustration  of  the  above-mentioned  vortex  evolntions  is  given 
in  Figure  5.7,  where  instantaneous  vorticity  isolines  are  shown  for  eight  character¬ 
istic  instances  dnring  a  flapping  cycle.  For  Uf/un  =  1/3  and  ojf/ujn  =  1/4,  it  is 
clear  that  the  enhanced  ESV  vortices  produce  an  oblique-shaped  TEV  vortex.  For 
ojf/ojn  =  1/4,  an  excessive  negative  camber  is  prodnced  at  stroke  reversal,  which 
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then  generates  a  high  suction  zone  on  the  lower  side  of  the  wing  leading  to  the  neg¬ 
ative  peak  in  the  CL  curve  seen  in  Figure  5.4.  The  Reynolds  number  effects  on  the 
temporal  evolution  of  the  lift  and  drag  forces  seen  in  Figure  5.4  can  also  be  observed 
in  Figure  5.8,  where  the  instantaneous  vorticity  isolines  are  shown  for  Re  =  250  for 
all  frequency  ratios.  Clearly,  as  the  viscous  damping  is  decreased,  the  system  dy¬ 
namics  system  ceases  to  be  periodic  and  the  important  vortices  are  stronger  and  are 
not  dissipated  as  quickly  as  seen  in  the  Re  =  75  case  (see  Figure  5.7).  For  the  case 
of  a  rigid  wing,  for  example,  the  LEV  from  a  given  stroke  interacts  with  the  shear 
layer  being  generated  in  the  next  stroke,  and  this  induces  a  premature  formation  of 
the  new  LEV.  This  process  is  not  periodic,  which  is  also  reflected  in  the  evolution  of 
lift  and  drag  forces.  A  similar  interaction  is  observed  at  ujjiujn  =  1/2  (see  last  three 
frames  in  Figure  5.8b).  Also,  in  all  flexibility  the  cases,  especially  for  uif/uin  =  1/3, 
strong  interaction  between  the  prohle  and  previous  ESV  is  seen  during  the  transla¬ 
tion  phase.  Regardless,  the  usual  components  of  the  flow  held,  namely  TEV,  LEV 
and  ESV  vortices  are  noticed  for  the  different  hexible  prohles  studied. 

5.1.4  Remarks  on  flexibility  effects 

Insect  wings  are  hexible  structures  that  undergo  large  displacements  and  de¬ 
formations  during  happing,  as  the  wing  structures  interact  with  the  surrounding 
how.  There  has  been  speculation  that  many  insects  hap  their  wings  at  frequencies 
close  to  the  natural  frequency  of  the  structure.  For  example,  analysis  of  Manduca 
sexta  wings  (see  for  example  references  [23],  [111])  has  shown  that  the  wings  hrst 
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Figure  5.7:  Instantaneous  vorticity  contours  over  one  period  at  Re  =  75.  Contours 
range  from  -10  (blue)  to  10  (red)  with  80  isolines.  Column  (a),  rigid  wing;  columns 
(b-d),  flexible  profile  with  ujf/ujn  =  1/2, 1/3,  and  1/4,  respectively. 
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Figure  5.8:  Instantaneous  vorticity  contours  over  one  period  at  Re  =  75.  Contours 
range  from  -10  (blue)  to  10  (red)  with  80  isolines.  Column  (a),  rigid  wing;  columns 
(b-d),  flexible  profile  with  uif  /ui^  =  1/2, 1/3,  and  1/4,  respectively. 
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natural  frequency  is  close  to  the  driving  frequency  in  normal  flapping  motion.  This 
suggests  that  insects  may  be  taking  advantage  of  a  structural  resonance  to  reduce  en¬ 
ergy  consumption  and  enhance  aerodynamic  performance.  Despite  the  signihcance 
of  such  a  hypothesis,  only  a  limited  number  of  studies  have  addressed  the  problem 
due  to  the  exceedingly  complex  fluid-structure  interactions  that  are  encountered  in 
experimental  or  numerical  work.  Although  three-dimensional  computations  of  the 
NavierStokes  system  coupled  with  a  wing  structural  system  are  within  the  reach  of 
todays  computers,  one  still  needs  to  develop  appropriate  mathematical  models  and 
tools  to  capture  all  important  phenomena  in  this  complex  system.  In  this  regard, 
the  present  study  extends  computational  work  that  has  been  conducted  before  with 
simplihed  two-dimensional  rigid  wings  to  include  the  effects  of  flexibility.  The  wing 
is  represented  by  two  rigid  links,  which  are  joined  at  the  center  by  a  pin  that  con¬ 
tains  a  torsion  spring.  The  kinematics  of  one  of  the  links  is  prescribed,  while  the 
motion  of  the  other  link  is  determined  by  the  fluidstructure  interactions.  Although 
a  fairly  wide  range  of  Reynolds  numbers  and  frequency  ratios  was  examined,  we 
found  that  the  computations  would  fail  at  forcing  frequencies  close  to  the  linear 
resonance  frequency  due  to  an  excessive  degree  of  flexibility.  This  limitation  is  due 
to  the  concentration  of  flexibility  at  a  discrete  point,  and  replacement  of  rigid  links 
with  elastic  links  modeled  as  elastic  beams  is  expected  to  help  in  overcoming  this 
limitation. 

As  mentioned  earlier,  the  two-link  structural  system  can  be  perceived  as  a 
double  pendulum  with  a  common  hinge.  In  particular,  when  one  of  the  link  motions 
is  prescribed,  the  other  link  behaves  as  a  pendulum  subjected  to  a  constraint  arising 
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from  the  prescribed  motion  and  complex  fluidstructure  interactions.  Equation  (5.2) 
does  resemble  the  equation  of  a  pendulum  driven  in  a  fluid.  A  straightforward 
perturbation  analysis  [71]  shows  that  the  structural  system  can  exhibit  non-linear 
resonances  at  cu/  =  l/3cc;„  and  ujf  =  3uJn-  These  resonances  are  expected  in  systems 
with  cubic  non-linearities;  for  example,  in  the  equations  governing  local  oscillations 
of  a  pendulum  about  an  equilibrium  position  and  elastic  systems  such  as  beams  [5] . 
Operation  of  the  flexible  wing  at  the  non-linear  superharmonic  resonance  Uf  = 
l/3ujn  was  found  to  be  benehcial  for  aerodynamic  performance.  Inclusion  of  fluid 
effects  will  give  rise  to  quadratic  non-linearities  and  additional  non-linear  resonances. 

For  the  specihc  set  of  kinematics  that  we  considered,  most  of  the  benehts  of 
having  a  flexible  wing  are  associated  with  the  stroke  reversal  phase  of  the  cycle. 
Especially  for  the  optimal  flexibility  cases  {ojf  =  l/Scu^),  the  strength  and  timing 
of  the  ESV,  as  well  as  the  dynamical  changes  of  the  wings  camber  due  to  struc¬ 
tural  deformations,  are  responsible  for  the  performance  enhancement.  The  overall 
enhancement  mechanism  is  analogous  to  that  produced  by  orientation  advancement 
in  rigid  wings  (see,  for  example,  reference  [107]).  It  is  noted  that  the  present  com¬ 
putations  cover  a  wide  range  of  frequency  ratios  and,  consequently,  wing  deflections 
range  from  a  few  degrees  to  very  large  values.  For  example,  in  the  case  of  highly 
stiff  wings  (e.g.  ujflujn  =  1/6),  the  maximum  deflection  angle  between  the  links  was 
about  11,  13  and  16  deg.  for  Re=75,  250  and  1000,  respectively,  while  for  highly 
flexible  wings  (e.g.  ujfjujn  =  1/2)  the  corresponding  numbers  were  67,  68  and  91 
deg.,  respectively.  In  insects,  the  wing  deformation  magnitudes  increase  as  the  body 
size  and  mass  increase,  and  it  is  conceivable  that  deformations  seen  in  this  study  at 
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the  aerodynamically  preferred  frequency  ratio  of  uif/uin  =  1/3  could  be  possible  in 
some  species.  On  the  other  hand,  for  small  insects  such  as  Drosophila  only  small 
magnitude  wing  deformations  have  been  observed.  The  computations  of  this  study 
show  that  as  the  wing  is  made  stiffer,  the  performance  enhancements  are  marginal 
when  compared  with  a  rigid  foil.  For  example,  at  Re  =  75  and  the  highly  stiff  case 
ojf/ojn  =  1/6,  Cl/Cd  is  approximately  6%  higher  than  that  of  a  rigid  foil.  For  the 
higher  Reynolds  numbers  Re  =  250  and  1000  there  is  actually  no  enhancement, 
and  the  performance  is  worse  than  that  obtained  with  a  rigid  foil.  The  above  re¬ 
sults  indicate  that  low  Reynolds  number  regimes  might  beneht  performance  even  at 
small  chordwise  distortions.  The  force  histories,  in  particular  for  the  low  Reynolds 
numbers,  appear  to  reach  a  periodic  steady  state  after  the  initial  transients  for  all 
of  the  frequency  ratios  that  were  considered,  suggesting  that  quasi-steady  models 
might  be  able  to  reproduce  this  behavior.  Such  models  have  been  reported  in  the 
literature,  and  have  been  adapted  for  flapping  flight  based  on  models  developed  for 
high  Reynolds  number  hxed  wing  aeroelasticity  studies  by  including  wing  rotation 
along  with  translation  ([32];  [34]),  and  forces  due  to  added  mass  ([89]).  In  a  more  re¬ 
cent  study  by  Wang  and  colleagues  [107]  the  unsteady  forces  from  experiments  and 
with  two-dimensional  computations  were  compared  with  the  quasi-steady  model 
predictions.  They  pointed  out  that  the  force  predictions,  which  were  made  by  using 
models  based  on  potential  flow  theory  [70]  for  a  constant  pitching  amplitude  and 
constant  translating  speed  wing,  deviated  substantially  from  the  experimentally  de¬ 
termined  unsteady  forces.  The  force  predictions  from  a  semi-empirical  model  based 
on  numerical  results  from  steady  translating  wings  at  a  fixed  angle  of  attack  were  in 
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broad  qualitative  agreement  with  the  unsteady  forces.  However,  detailed  compar¬ 
isons  revealed  that,  depending  on  the  kinematics,  the  unsteady  effects  can  reduce 
the  total  lift  by  a  factor  of  2  to  3.  In  the  present  case,  due  to  the  wings  flexibility, 
the  identihcation  of  the  quasi-steady  contributions  is  more  complex  as  additional 
new  states  have  been  included. 

5.2  An  Insect  at  free  longitudinal  flight 

While  aerodynamics  of  insects  and  other  flying  animals  has  received  consid¬ 
erable  attention  over  the  past  decade,  the  same  is  not  true  of  response  stability  of 
flying  animals.  Currently  there  is  a  lack  of  quantitative  understanding  of  the  sta¬ 
bility  of  flying  insects.  In  this  section  we  will  present  preliminary  results  from  our 
effort  towards  the  development  of  a  comprehensive  framework  that  will  enable  us  to 
conduct  sophisticated  numerical  experiments  to  understand  and  assess  the  response 
stability  of  insect-like  systems  in  different  disturbance  environments. 

5.2.1  Multi-body  insect  model  kinematics 

The  insect  model  is  composed  of  one  pair  of  rigid  wings,  RW  and  LW ,  hinged 
to  a  rigid  body,  RBI,  which  respresents  the  insect’s  head  and  thorax,  and  another 
body,  RB2,  representing  the  insect’s  abdomen.  The  latter  is  also  articulated  to 
RBI  (See  Figure  5.9).  The  coordinate  transformations  among  the  different  reference 
frames  and  with  respect  to  the  inertial  frame  N  are  dehned  in  terms  of  Euler  angle 
sequences.  The  kinematic  analysis  for  the  given  system  starts  by  selecting  the 
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RB1 


Figure  5.9:  Four  rigid  body  model  for  the  full  insect.  Rigid  Body  RBI,  head-thorax, 
is  represented  by  the  tracking  frame  B  attached  to  its  center  of  mass.  Rigid  body 
RB2,  the  abdomen,  is  hinged  to  RBI  at  point  D  where  its  respective  body  frame 
is  defined.  The  wing  bodies  RW  and  LW  are  articulated  to  RBI  at  points  R  and 
L  respectively,  where  their  tracking  frames  are  defined. 

sequences  that  will  define  the  transformations  among  the  reference  systems  N,  B,  D, 
R  and  L.  The  transformation  between  the  inertial  system  N  and  the  RBI  tracking 
system,  B,  defined  by  the  unitary  vectors  bi,  b2,  bs,  is  given  by  a  180°  rotation 
with  respect  to  the  axis  given  by  rii,  followed  by  a  3-2-1  Euler  angle  sequence  with 
angular  coordinates  xjj  (yaw  respect  to  63),  6  (pitch  respect  to  62)  and  0  (roll  respect 
to  e”),  as  shown  in  Figure  5.10.  This  is  the  typical  flight  mechanics  coordinate 
transformation  and  yields  a  singularity  for  pitch  angles  multiple  of  ti/2. 
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Figure  5.10:  Transformation  from  inertial  reference  frame  N  and  RBI  body  frame 
B.  A  180°  rotation  respect  to  fii  is  given,  followed  by  a  3-2-1  Euler  angle  sequence 
with  angular  coordinates  xjj  (yaw  respect  to  es),  6  (pitch  respect  to  62)  and  0  (roll 
respect  to  e'() 


The  transformation  matrix  from  N  io  B  is 


Tsat]  — 


ddcijj 


—  c9sxjj 


s9 


s(j)s6cijj  —  ccl)S'ip  —s(j)s6s'ijj  —  c(j)cip  —s(j)c6 


ccpsOcijj  -|-  s(j)S'ijj  —C(j)s6s'ijj  -|-  s(j)cxjj  —c(j)c6 


(5.9) 


where  s()  and  c()  are  the  sin()  and  cos()  functions  respectively.  The  transformation 
from  system  B  to  the  system  D  local  to  the  abdomen  is  given  by  two  rotation  angles: 
angle  ip2  with  respect  to  ba,  and  62  with  respect  to  the  intermediate  axis  b2.  The 
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corresponding  transformation  matrix  is: 

06^2  C02  C02#2  -SO2 

[Tdb]  =  —  5-02  C'02  0  (5.10) 

s92Cll)2  56*2  #2  062 

Next  the  transformation  between  the  system  B  and  the  system  R  local  to  the  right 
wing  is  done  in  the  following  steps: 

1.  90°  rotation  with  respect  to  bs,  followed  by  a  90°  rotation  with  respect  to 

"  // 

2.  /dg  degrees  rotation  with  respect  to  bg 

-  /// 

3.  An  3-2-1  Euler  angle  sequence  of  angles  '^3  with  respect  to  bg  (beating  angle 
in  the  stroke  plane,  dehned  here  relative  to  the  body  coordinate  system),  9^ 

^  IV  -  V 

with  respect  to  b2  (out  of  stroke  plane  angle),  and  03  with  respect  to  bg 
(angle  that  modihes  the  angle  of  attack  of  the  wing). 

The  resulting  transformation  matrix,  [T^^],  is  dehned  as  follows: 

cflgsbasdg  -  s6»3Cd3  C03C03  cOaSTpsCfi'^  +  503  503 

{a4s93  +  ai)sf3'^  +  S(j)3c93cP'^  a2S03  —  ag  (q;4S03 -I- Q;i)c/33  —  (5-11) 

(03  S03  -  02)5/33  +  01503-1-04  {^3393  -  a2)  -  c4>3c93Sl3'3 

where  the  R  system  is  dehned  by  the  versor  triad  ri,r2,r3  (see  Figure  5.11  for 
reference  on  the  diherent  coordinate  systems).  Here,  ai  =  0:2  = 

0:3  =  and  0:4  =  The  transformation  from  frame  B  to  frame  L  is  very 

similar  to  the  one  described  above: 

1.  —90°  rotation  with  respect  to  b3,  followed  by  a  —90°  rotation  to  respect  to  bg 
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2.  degrees  of  rotation  with  respect  to  b^.  Assuming  that  the  pitch  angle  of 
body  B  is  zero,  the  angle  (5'^  is  related  to  the  stroke  plane  angle,  (3,  reported 
in  many  studies  ([108],  [7])  by  the  equation  =  90°  —  (3. 

^  /// 

3.  Euler  angle  sequence  3-2-1  of  angles  1^4^  with  respect  to  bg  ,  64  with  respect  to 

-IV  -  V 

b2  ,  and  04  with  respect  to  b^ . 

In  the  same  manner  the  transformation  matrix  (04,  ^4, 04;  04)]  can  be  written 
as: 

C04S04S/34  —  S04C/34  —  C04C04  —  C04S04C/34  —  5^4  504 

(7456*4 -I- 7i) 504 -I-  504  004  004  ”72504  “t  73  - (74 504  “t  7l ) 004 -|-  504 004 504  (5-12) 

(73504  -  72)504  -1-  004004004  ”71504  -  74  ”(73504  -  72)  004 -|-  004  004504 

The  L  system  is  dehned  by  the  versor  triad  li,  12,  Is,  and  71  =  c(f)4C'ip4,  72  =  s04C04, 
73  =  c(j)4Silj4,  74  =  S04S04.  The  value  that  must  be  given  to  the  angular  variable  03 
on  the  right  wing,  in  order  to  conserve  the  same  stroke  plane  as  the  one  present  for 
the  left  wing  is  03  =  —04. 

In  order  to  produce  symmetric  flapping,  the  relation  among  the  right  wing 
Euler  angle  sequence  and  left  wing  angular  variables  must  be: 

03 (t)  =  -04 (t)  ;  6^3 (t)  =  04 (t)  ;  03  (t)  =  -04(0  (5-13) 

Once  the  transformations  among  coordinate  systems  and  the  positions  of  the 
origins  of  systems  D,  R  and  L  are  given  with  respect  to  the  origin  of  the  thorax, 
B,  the  positions  of  all  points  in  the  different  component  bodies  can  be  computed  in 
terms  of  the  inertial  frame  N.  Since  the  rotation  matrices  form  an  orthonormal  set 
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Figure  5.11:  Reference  frames  for  rigid  body  system, 
for  the  points  of  body  B  we  have: 

<  =  R"  +  |Ta,b|  r«/B  (5.14) 

where  =  [x(t)  y(t)  z(t)],  [Ttvb]  —  [T^tv]^,  represents  the  position  of  a  generic 
point  in  body  B  with  respect  to  the  origin  of  the  inertial  frame  in  terms  of  the  hi, 
h2,  ha  triad,  and  ^p/B  position  vector  from  the  origin  of  frame  B  to  point  p  in 

local  bi,  b2,  ba  coordinates.  A  similar  expression  can  be  written  for  the  remaining 
bodies: 


+  [Tats]  +  [Tatd]  {point  in  the  abdomen)  (5.15) 

=  Rs  +  [Tats]  r^/B  +  [Tair]  {point  in  right  wing)  (5.16) 

^PL  =  +  [TNB\rL/B  +  \^NL\rp^/L  {point  in  left  wing)  (5.17) 
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To  complete  the  kinematic  analysis  of  this  mechanical  system  the  angular  velocities 
and  accelerations  of  the  different  bodies  in  terms  of  the  time  derivatives  of  the  Euler 
angles  must  be  obtained.  In  hgure  5.10,  the  angular  velocity  of  frame  B  with  respect 
to  the  inertial  reference  frame  N  can  be  expressed  in  terms  of  the  time  derivatives 
of  the  Euler  angles: 

=  Tpes  +  9e2  +  <j)e”  (5.18) 

where  the  ea,  62,  e”  unit  vectors  in  terms  of  the  versors  of  the  B  frame  are: 


■ 

■ 

<  > 

1 

0 

-sO 

0 

N,  ,B  = 

Ll>b 

0 

ce 

c6s(j) 

< 

9 

0 

—  S(f) 

c9c(f) 

0 

V.  / 

es  =  — s6*bi  +  c6s(l)h2  +  cOccjyhs  ;  62  =  c0b2  —  s^bs  ;  e”  =  bi  (5.19) 
Substituting  ((5.19))  into  ((5.18))  and  rearranging  into  matrix-vector  form  we  get: 


(5.20) 


which  is  the  angular  velocity  of  frame  B  respect  to  N  in  terms  of  the  unit  vectors 
bi,  b2,  ba-  We  denote  the  matrix  in  the  above  equation  as  and  0^  = 

[0  6  0]^.  Considering  the  vector  identity  x  =  0,  the  corresponding 

angular  acceleration  is  given  by: 

''af  =  ^  +  [«B  (5.21) 

The  velocity  and  acceleration  of  a  point  p  on  RBI,  in  terms  of  the  unit  vectors  of 
frame  B  are: 


B  ^  d^ 


N 


d 


—  {Rb)b+  ^B^^p/B 


(5.22) 


dB 


:Rb)b  +  X  r^/B  +  X-  u.-  X  r^/B  (5-23) 
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The  angular  velocity  of  frame  D  with  respect  to  frame  B  is  +  ^2102) 

and  in  terms  of  the  B  frame  decomposition: 

-#2  0 
C1p2  0 
0  1 

Then,  the  angular  velocity  and  angular  acceleration  of  frame  D  with  respect  to  the 
inertial  frame  N  is: 


92 

^2 


(5.24) 


B,  ,D  = 


'Ll) 


B 


(5.25) 

N^D  ^  ^  rN^B  ^  ^  ^  ^  ^  /B^BN  ^  N ^B  ^ 

at  ^  '  at  ^  '  at  ^  ' 

Then,  using  equation  (5.24),  ^  =  [-^B +  the  absolute 

angular  acceleration  of  frame  D  in  terms  of  the  B  frame  unit  vectors  is: 


N „D  N  „B 


B„D 


N,  .B 


=  a%+  x 


B,  .D 


'■B 


B 


'B 


(jJ 


B 


(5.27) 


The  linear  velocity  and  acceleration  of  a  point  on  RB2  can  be  found  as  follows: 


PD 


I  N  .B  ^  B 

b)b+  X  r 


D/B 


+  X 


X  r,^ 


PD 


Pd/D 


+  (jJ^  X 


PdID 

,B 

Pd/D 


(5.28) 

(5.29) 


where  a^  is  the  acceleration  of  the  D  frame  origin,  and  is  computed  from  equa¬ 
tion  (5.23)  for  this  point  in  RBI. 

For  the  right  wing,  the  angular  velocity  dehnition,  and  its  expression  in  terms 
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of  frame  B  are: 


•  -///  •  -/y  •  -y 

'j/’sbg  +  631)2  +  03bi  ,  and 


0638^)33(33  -  3630(33 

0^)38(33 

Cp3 

<  > 

03 

B,  ,R  _ 

U)b  — 

0630^)3 

-#3 

0 

< 

^3 

0638^)30(33  -f-  3633(33 

0^)30(33 

-s(33 

03 

^  y 

(5.30) 


(5.31) 


respectively.  Then  the  corresponding  absolnte  angnlar  velocity  and  accelera¬ 
tion  can  be  readily  obtained  from  expressions  analogous  to  equations  (5.25)- 
(5.27).  Then  the  linear  velocity  and  acceleration  for  a  point  pn  laying  on  the  right 
wing,  RW,  are: 


N 


< 


d 


B 


=  ^  (Rb)b  +  X 


-h  X 


=  -I-  X 


vr/R 


+  (jj^  X 


pr/R 

.B 

pr/R 


(5.32) 

(5.33) 


Equations  similar  to  (5.25)-(5.27)  and  (5.32)-(5.33)  can  de  derived  for  the  kinematics 
of  the  left  wing  in  terms  of  frame  B.  The  angular  velocity,  for  example,  relative  to 
B  is: 


or  •  -///  ■  ^  IV  .  -y 

=  '^4b3  -|-  6*4b2  -|- 


c6*4  504504  —  56*4  C04 

C04504 

C04 

<  > 

04 

B,  ,L  _ 

Ll>b  - 

—  0640(4 

504 

0 

< 

04 

—  {06481(40(34  -|-  56^4504) 

—  C04C04 

504 

04 

<  y 

(5.34) 


(5.35) 


5.2.2  Kinematics  of  longitudinal  flight 


In  the  present  study  we  will  limit  ourselves  to  longitudinal  flight,  where  the 
insect’s  motion  is  simplified  to  longitudinal  mechanics.  Such  a  flight  mode  is  possible 
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when:  i)  a  symmetric  wing  motion  is  used  and  the  lateral  force  and  moments  are 
very  small  compared  to  their  longitudinal  counterparts;  ii)  the  mass  distribution  is 
symmetric  and  has  the  same  symmetry  plane  as  the  model  geometry.  In  such  case 
the  the  coordinates  for  RBI  reduce  to  the  displacements  x{t)  and  z{t)  along  the 
symmetry  plane,  fii  —  113,  and  the  6{t)  angle  with  respect  to  the  body  axis,  b2, 
which  is  aligned  at  all  times  with  the  — n2  direction.  The  transformation  matrix 
from  frame  N  to  B  and  the  angular  velocity  vector  reduce  to: 

cos{6)  0  sin{6) 

[Tbv]=  0-10,  =  (5.36) 

sin{6)  0  —cos{6) 

The  coordinates  that  describe  the  relative  orientation  of  RB2  with  respect  to  RBI 
reduce  to  the  variable  62.  The  corresponding  tansformation  matrix  from  frame  B 
to  D,  and  angular  velocity  are: 

003(62)  0  —sin(62) 

[Tdb]=  0  10,  =  (5.37) 

sin(62)  0  003(62) 

The  rest  of  the  kinematic  analysis  follows  the  same  path  as  in  the  previous  section. 

5.2.3  Dynamics  of  longitudinal  flight 

Here  we  derive  the  equations  of  motion  for  the  longitudinal  dynamics  using  the 
classical  Lagrange’s  equations  approach  [14].  We  dehne  the  generalized  coordinates 
x(t),  z(t),  6(t)  and  6*2 (t)  to  be  force  driven,  and  all  angular  variables  that  describe 
the  relative  motion  of  the  wing  with  respect  to  RBI  to  be  prescribed  in  time.  We 
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first  need  to  define  the  mass  inertia  properties  associated  with  each  of  the  four  rigid 
bodies.  For  RBI  they  are:  and  Isyy^  the  total  mass  and  mass  moment  of  inertia 

with  respect  to  the  body  axis  given  by  b2  (see  hgure  5.12a).  We  make  the  additional 
assumption  that  the  origin  of  the  B  frame  is  at  the  same  location  as  the  center  of 
mass,  CMrbii  of  body  RBI.  All  mass  moments  of  inertia  are  dehned  with  respect 
to  the  corresponding  body  center  of  mass  in  local  coordinates.  For  RB2\  tud  and 
iDyyt  are  the  total  mass  and  mass  moment  of  inertia  with  respect  to  the  body  axis 
passing  through  CMbb2  and  is  parallel  to  d2,  which  is  in  longitudinal  motion  the 
same  as  b2.  The  center  of  mass  of  RB2  is  located  on  the  symmetry  plane,  di  —  da, 
at  positions  xqd  and  zcd  from  the  origin  of  frame  D  (see  hgure  5.12b). 

The  right  wing,  RW,  is  on  f i  —  fs  plane  and  in  hnite  thickness  cases  the  wing  is 
assumed  to  be  symmetric  with  respect  to  this  plane.  The  required  inertia  properties 
in  this  case  are:  uir,  Irxxi  iRyy,  Irzzi  Irxz,  the  wing  mass  and  moments  of  inertia 
with  respect  to  axes  passing  through  CMrw  and  parallel  to  fi,  f2  and  fa,  and  the 
product  of  inertia  on  the  fi  —  fa  plane  respectively  (see  hgure  5.12c).  The  center  of 
mass  of  RW  is  located  at  the  local  coordinates  xcr  and  zcr  from  the  origin  of  the 
R  system.  The  left  wing,  LW,  has  similar  properties  as  the  right  wing:  m^,  Ilxxi 
^Lyyi  Ilzzi  Irxz-  Due  to  the  symmetry  assumption  the  inertia  properties  of  the  left 
wing  are  equal  to  the  ones  on  the  right  wing.  The  position  of  the  center  of  mass  of 
LW  is  given  by  the  xcL)  zcl  local  coordinates. 

Next  we  will  dehne  the  kinetic  and  potential  energies  for  each  rigid  body  in 
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Figure  5.12:  Inertia  properties  and  reference  frames  for  bodies:  (a)  RBI,  (b)  RB2, 
and  (c)  RW.  Properties  of  the  left  wing  LW  are  analogous  to  RW.  The  degrees  of 
freedom  of  the  system  are  the  variables  x{t),  z{t),  6{t)  and  6*2(t). 

the  model.  In  particular  the  kinetic  energy  can  be  dehned  as: 

Tb,  =  2  rriBi^cMg.  +  ^  [Ib,]  ^ (5.38) 

where  Bi  is  any  of  the  RBI,  RB2,  RW  or  LW  bodies;  ^cmb-  f^e  absolute 
velocity  vector  of  the  center  of  mass  of  Bi,  is  the  absolute  angular  velocity  of 

Bi]  and  ms.,  [J^j  are  the  mass  and  mass  inertia  matrix  of  Bi.  These  vectors,  as 
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we  discussed  in  the  previous  section,  can  be  expressed  in  terms  of  the  B  frame  unit 
vectors  bib2b3.  The  inertia  matrix  for  each  body  also  needs  to  be  expressed  in  terms 
of  these  versors.  The  potential  energy  for  each  body  is  given  by  the  gravitational 
held: 

Vsi  =  9  zcMBi  (5.39) 

where  g  is  the  gravitational  acceleration  and  zcmb-  is  the  vertical  position  of  the 
center  of  mass  coordinate  along  the  Newtonian  unit  vector  113  for  body  Bi. 

In  the  simplest  case  of  RBI  the  resulting  kinetic  and  potential  energy  functions 

are: 

TRBi  =  ]^rnB{x{tY  +  z{tf) +  ]^lByyO{tf,  Vrbi  =  niB  g  z{t)  (5.40) 

Then  applying  Lagrange’s  equations  on  the  Lagrangian  function  Lrbi  =  Trbi  — 
Vrbi,  using  x(t),  z(t),  6{t)  and  6'2(t),  as  generalized  coordinates  we  get  the  following 
system: 
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where  Qbx,  Qbz,  Qbo  and  QBe2  are  the  generalized  huid  forces  acting  on  the  gen¬ 
eralized  coordinates.  We  see  that  the  62  equation  has  no  contributions  from  RBI 
(Qb02  is  zero).  The  contributions  to  the  62  equation  come  from  RB2.  Note  that 
for  RB2  there  is  a  component  of  the  potential  energy  arising  from  a  linear  torsional 
spring  Kr  located  at  point  D,  which  exerts  a  restituting  moment  proportional  to  62 
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(figure  5.12b).  A  torsional  damper  could  also  be  added.  The  computation  of  fluid 
forces  and  moments  for  each  of  the  immersed  bodies  in  the  model  is  done  using  the 
approach  in  Chapter  2,  which  enables  us  to  compute  the  force  =  [Fbx,  Fsy,  Fbz]^ 
and  moment  with  respect  to  the  center  of  mass,  M;^  =  [Mbx,  MBy,  Mbz^ ■,  of  RBI 
in  the  N  frame.  In  terms  of  the  B  frame  these  forces  are  simply  =  [Tbn]Fb 
and  =  [Tbn]  Mb-  Then,  the  generalized  fluid  forces  are  computed  from  [14] 


Qbx,  =  F^ 


^^RBl 


+  M^ 


J  =  l,2, 


(5.42) 


where  Xj  represents  any  of  x,  ;2,  6*,  62-  The  resulting  generalized  fluid  forces  on 
RBI  are:  Qbx  =  Fbx,  Qbz  =  Fbz,  Qbo  =  -MBy,  Qb02  =  0-  The  procedure  for 
the  equations  for  RB2,  RW  and  LW  is  similar,  where  the  corresponding  velocities 
of  the  centers  of  mass  and  absolute  angular  velocities  for  each  of  these  bodies  is 
utilized.  The  problem  becomes  intractable  to  derive  analytically,  specially  for  the 
wings  which  undergo  three-dimensional  motion  given  by  the  beating  and  ont  of 
stroke  plane  angles  and  6^3,  6^4  as  well  as  the  angles  of  attack  03,04.  For  this 

reason  the  derivation  of  the  eqnations  for  each  body  was  done  using  the  symbolic 
manipulation  software  Mathematica©. 


Testing  the  consistency  of  the  model 

To  verify  the  above  model  two  test  cases  were  setnp,  where  the  effects  of  the 
flnid  flow  are  neglected.  In  the  hrst  problem  the  system  is  composed  by  RBI  and 
RB2  only.  The  degrees  of  freedom  in  this  case  are  the  fonr:  x,  z,  6  and  62-  This 
system  undergoes  two-dimensional  motion  and  is  equivalent  to  the  two-link  model 
presented  in  Section  5.1.1.  If  the  inertial  properties  and  initial  conditions  are  selected 
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to  be  the  same,  then  both  sets  of  equations  should  provide  the  same  response.  A 
comparison  between  the  two  models  is  shown  in  Figure  5.13  for  the  following  initial 
conditions:  for  the  RB1-RB2  computation,  9o  =  —0.5,  620  =  1,  =  xdbcos{6o) 

and  Zo  =  0;  for  two-link  model,  60  =  —0.5  -|-  tt,  620  =  1,  a;o  =  0  and  Zo  =  rii,sin{9o). 
As  it  can  be  seen  the  evolution  of  the  resulting  energies  in  time  for  the  both  models 
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Figure  5.13:  Energy  plots  as  a  function  of  integration  time  for  comparison  among 
RB1-RB2  and  two-link  model:  (a)  (blue)  kinetic  energy  of  RBI,  (red)  kinetic  energy 
of  RB2,  (•)  corresponding  kinetic  energies  for  plates  B  and  A  of  the  two-link  model; 
(b)  (red)  total  energy  of  the  system  RB1-RB2,  (•)  total  energy  of  the  two-link 
model. 

match  perfectly.  The  maximum  difference  for  an  integration  of  15  computational 
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units  (letot)  was  of  the  order  of  the  tolerance  imposed  on  the  numerical  integration 
(about  10“^°).  Excellent  agreement  was  also  found  on  the  horizontal  and  vertical 
positions  of  the  centers  of  mass  of  each  body  (hgure  5.14).  These  results  conhrm 
that  the  contributions  of  RBI  and  RB2  to  the  four  rigid  body  system  in  longitudinal 
dynamics  are  correctly  computed. 


t 


t 


Figure  5.14:  Time  variation  of  positions  of  centers  of  mass:  (a)  Horizontal  positions, 
(blue)  RBI,  (red)  RB2,  (•)  corresponding  positions  for  plates  B  and  A  of  the  two- 
link  model;  (b)  Vertical  positions,  bodies  are  identihed  as  in  (a). 

The  second  test  was  performed  on  the  four  rigid  body  model,  using  harmonic 
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kinematics  on  the  wing’s  stroke  angle  and  angle  of  attack.  The  oscillation  amplitude, 
forcing  frequency,  and  also  the  inertia  and  stiffness  properties,  were  the  ones  outlined 
in  the  following  sections.  The  equations  of  motion  were  integrated  using  a  fourth 
order  Runge-Kutta  method  in  Matlab©  (function  ode45).  A  zero  gravity  force  was 
prescribed.  As  there  are  no  external  forces  applied  to  the  system  throughout  the 
calculation  the  position,  xcmt^UcMtt  ^cMt^  of  fhe  center  of  mass  of  the  four  rigid 
body  system: 


XcMt 

UCMt 

ZcMt 


rriBXcMjiBi  +  f^DXCMjiB2  +  ^R^CMbw  +  ^LXcMlw 
rriB  +  rriD  +  niR  +  mi 
RRrVCMbw  +  ^LVCMlw 
rriR  +  rriL 

RRbZcMrbi  +  ^dZcMrb2  +  ^R^CMrw  +  ^RZcMlw 
rriB  +  m-D  +  RRr  +  RRl 


(5.43) 


is  expected  to  be  stationary  in  the  Newtonian  frame  axes,  hi,  112,113.  The  time 
integration  was  performed  over  20  flapping  cycles  and  the  integrator  absolute  er¬ 
ror  tolerance  was  set  to  10“^^.  The  maximum  variations  of  the  xcmt^UcMt^  ^cmt 
variables  in  time  were  found  to  be  less  than  2.5e“^^,  which  is  of  the  order  of  the 
tolerance  level  imposed  on  the  numerical  integration.  Different  values  of  the  integra¬ 
tion  tolerance  provided  corresponding  variation  values  for  xcmtjVcmt^  ^cmt-  This 
test  showed  that,  although  the  generalized  coordinates  x(t),  z(t),9(t),92(t)  of  the 
model  did  undergo  important  variation,  there  is  no  spurious  momentum  introduced 
to  the  system. 
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5.2.4  Hovering  of  a  Musca  Domestica  model  at  Re  =  500 

In  this  section  we  will  present  simulations  for  a  fly  model  hovering  at  Reynolds 
number  of  500.  Both  prescribed  kinematics  and  FSI  of  free  flight  are  considered.  The 
model  was  created  using  CAD  software  from  Musca  Domestica  digital  images.  The 
different  bodies  can  be  seen  in  hgure  5.12.  The  reference  length  is  the  wingspan, 
b  (one  wing  hinge  to  tip  distance).  The  head-abdomen  length  of  the  model  is 
Lta  =  1.036.  The  wings  have  rounded  edges  and  a  thickness  of  0.0256.  A  simple 
set  of  wing  kinematics,  representative  of  Diptera  wing  motion  is  prescribed  in  all 
simulations.  For  the  right  wing,  for  example,  we  set: 

03(6)  =  A^sin  {ujft  -h  a) 

O^it)  =  0  (5.44) 

03(6)  =  Iprn  +  A^COS  {U ft)  , 

where  =  7r/4  is  the  amplitude  of  the  angle  of  attack,  and  a  =  71/6  (advanced 
rotation)  is  the  phase  .  The  mean  stroke  angle  is  ipm  =  — 7r/36,  and  the  stroke 
amplitude  =  557r/180.  The  reference  velocity  Ur  used  is  the  mean  wing  tip 
velocity  given  by  Ur  =  '•p^.mean^i  where  03mean  =  ‘lAfjUif/Ti.  The  Reynolds  number 
is  Re  =  URh/v  =  500.  The  RBI  orientation  angle  on  the  prescribed  kinematics 
simulation  was  set  to  6^  =  7r/3.  Using  a  stroke  angle,  03  =  p6/6,  the  stroke  plane  was 
aligned  with  the  horizontal  plane.  The  values  given  to  the  angles  of  the  left  wing 
were  such  that  symmetric  flapping  occurs. 

Prescribed  kinematics  simulation 
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First,  a  calculation  of  a  tethered  fly  with  with  prescribed  kinematics  was  con¬ 
ducted.  The  relative  angle,  62,  between  RBI  and  RB2  was  set  to  zero.  The  model 
was  placed  in  the  center  of  a  cubic  domain  of  dimensions  7b  x  7b  x  76,  with  the 
origin  of  frame  i?  at  a  distance  z  =  4b  from  the  bottom  wall.  No-slip  boundary 
conditions  were  set  on  the  x  and  directions,  and  no  penetration  on  the  y  direction, 
where  the  Poisson  equation  was  assumed  to  have  a  periodic  solution.  Level  0  on 
the  AMR  grid  consists  of  8  x  8  x  8  blocks  with  16^  cells.  Three  levels  of  rehne- 
ment  were  used,  and  together  with  dynamic  adaptivity,  the  cell  size  around  all  solid 
boundaries  was  Ax  =  Ay  =  Az  =  0.00686.  With  this  resolution  approximately  150 
points  were  clustered  along  each  wing  span.  Mesh  adaptation  was  performed  every 
10  timesteps.  The  grid  was  refined/derehned  based  on  two  criteria:  the  presence 
of  a  solid  boundary  and/or  the  magnitude  of  the  vorticity  vector  on  a  block.  A 
block  was  rehned  when  the  maximum  vorticity  magnitude  within  it  was  larger  than 
5.5Uii/b,  and  a  set  of  children  where  derehned  when  it  was  less  than  4.817ij/6.  The 
problem  was  integrated  at  constant  CFL  number,  for  a  time  span  of  four  flapping 
cycles. 

In  order  to  examine  the  symmetry  of  the  flow  and  corresponding  loads  we  eval¬ 
uate  the  maximum  total  force  in  the  Eulerian  reference  frame.  The  corresponding 
force  coefficients  are  dehned  as: 

(~imax  ^  mao;  +  F^iRB2it)  +  F^jRwif)  +  F^^Lwif)) 

\PfUlS^ 

where  Xj,  i  =  1,  2,  3  represents  the  Eulerian  x,  y  and  2:  coordinates,  and  S^,  =  0.2746^ 
is  the  one-wing  planform  area.  The  forces  F^^lt)  are  the  hydrodynamic  forces  for 
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each  of  the  bodies  in  the  model. 


The  maximum  fluid  force  coefficients  in  the  Eulerian  frame  are  (7™“®  =  3.35, 
(jmax  _  5_3g-3^  (7™“^  =  2.13.  The  ffuid  force  on  the  y  direction  normal  to  the 

model  symmetry  plane  is  0.25%  the  value  of  the  vertical  force,  which  demonstrates 
the  validity  of  the  longitudinal  flight  assumption  for  this  particular  set  of  wing 
kinematics.  In  hgure  5.15  the  time  variation  of  Cx  for  each  of  the  four  rigid  bodies 
is  shown.  Both  RW  and  LW  results  are  identical,  and  bodies  RBI  and  RB2 
contribute  a  maximum  Cx  that  is  6%  of  the  amount  produced  by  the  wings.  This 


Figure  5.15:  Time  variation  of  Cx  coefficient  for  different  bodies  on  prescribed  sim¬ 
ulation:  (blue)  RBI,  (dashed  magenta)  RB2,  (red)  RW,  and  (•)  LW. 

result  is  important  since  it  points  to  the  need  of  modeling  the  thorax-head  and 
abdomen  bodies  in  some  high  resolution  simulations  of  hovering.  The  lift  and  drag 
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coefficients  for  this  normal  hovering  simulation  is  computed  for  RW  by 


CluW 

CduW 


FzRwjt) 

—sgn 


FiRw(t)  cos  +  F,Rw(t)  sin  (j’sd)) 

\P,UIS^ 


(6.46) 


(5.47) 


and  similar  expressions  can  be  used  for  RW .  Also  the  lift  coefficient  can  be  computed 
directly  for  RBI  and  RB2.  These  results  can  be  seen  in  hgure  5.16.  The  mean  lift 
and  drag  coefficients  computed  in  periods  3  and  4  for  the  right  wing  are  = 

0.42  and  Cd^w  ~  0.77.  Same  coefficients  were  found  for  LW .  The  maximum 


Figure  5.16:  Time  variation  of  Cl  and  Cl,  coefficients  for  different  bodies  on  pre¬ 
scribed  simulation:  (blue)  RBI,  (dashed  magenta)  RB2,  (red)  RW ,  and  (•)  LW . 

contribution  to  lift  provided  by  RBI  and  RB2  is  about  4%  of  the  lift  given  by  both 
wings. 
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Free  longitudinal  flight  simulation 

To  be  able  to  conduct  simulations  of  free  flapping  flight  we  hrst  need  to  get 
estimates  for  the  inertia  properties  for  each  of  the  rigid  bodies  in  our  model.  Since 
there  are  no  readily  available  data  in  the  literature,  we  present  here  an  approach 
to  obtain  realistic  estimates  for  the  model  we  use.  We  start  by  setting  the  ratio  of 
inertia  to  aerodynamic  force  we  expect  to  obtain  on  the  course  of  the  simulation.  An 
estimate  of  this  ratio  is  given  by  comparing  the  mean  inertia  force  due  to  tangential 
acceleration  of  the  wing  center  of  mass  Fin  =  ‘2^mwXcRA^u^/n  to  the  mean  drag 
force  Dw  from  the  prescribed  kinematics  calculation  above.  In  a  manner  similar  to 
the  2D  airfoil  calculations  presented  earlier,  we  assume  the  regime  is  one  of  high 
fluid-structure  coupling  and  Fin/D^  —  1.3.  The  inertia  effects  are  expected  to  be 
higher  due  to  the  contribution  of  the  mass  moments  of  inertia  in  three  dimensional 
wing  motion.  From  this  ratio,  the  mean  drag  force  coefficient  previously  obtained, 
and  the  wing  geometry  we  can  derive  the  value  of  rriw  =  0.2.  All  inertia  properties  are 
made  dimensionless  using  the  fluid  reference  variables  pf  and  b.  On  the  other  hand, 
from  the  wings  CAD  model  the  volume  contained  by  the  wing  can  be  computed. 
Assuming  a  constant  density  of  the  wing  Pw/Pf  —  30  the  mass  estimated  before 
can  be  obtained.  Then,  all  mass  inertia  properties  for  the  wings  can  be  computed 
performing  numerical  integration  along  their  volume 

Iwxx  =  (y"  +  2^")  Iwyy  =  +  Z^)  dQ  (5-48) 

Iwzz  =  {x^  +  y'^)dVt  Inixz  =  X  zdVt 

here  w  refers  to  any  of  RW  and  LW .  All  integrations  are  performed  in  the  local 
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Body  i 

tub. 

^Bixx 

^Biyy 

^Bizz 

^Bixz 

XcBi 

^CBi 

RBI 

2.323 

— 

6.646e“^ 

— 

— 

0.0 

0.0 

RB2 

1.547 

— 

2.448e-2 

— 

— 

-0.180 

-1.984e-2 

RW 

0.2 

1.741e-3 

1.326e-2 

1.154e-2 

3.962e-^ 

0.500 

-9.923e-2 

LW 

0.2 

1.741e-3 

1.326e-2 

1.154e-2 

3.962e-^ 

0.500 

-9.923e-2 

Table  5.2:  Dimensionless  inertia  properties  for  bodies  of  four  rigid  body  model, 
computed  in  local  coordinate  systems.  xcBi  and  zcBi  are  the  local  coordinates  of 
the  center  of  mass  of  body  5*  respect  to  the  origin  of  the  local  frame  of  reference. 


coordinate  system  centered  on  the  wings  center  of  mass  location.  The  numerical 
integration  was  performed  by  meshing  the  wings  volume  with  linear  tetrahedra  ele¬ 
ments,  and  performing  a  finite  element  sum  of  4  point  quadratures  on  the  elements. 
This  quadrature  and  the  linear  interpolation  functions  used  on  the  spatial  variables 
ensures  that  the  integrations  are  exact  on  the  discretized  volumes.  Grid  rehnement 
was  performed  to  ensure  that  the  inertia  properties  were  converged  up  to  10“®. 

To  dehne  the  mass  of  the  thorax-head  and  abdomen  we  assume  that  one  wing 
mass  in  our  model  is  5%  of  the  mass  of  RBI  and  RB2  together.  Wu  et.  al  [112] 
report  this  percentage  to  be  1%  for  Hoverflyes  and  6%  for  the  Manduca  Sexta 
hawkmoth.  Then,  the  density  of  RBI  and  RB2  was  found  to  be  Pb/ Pf  —  50.  The 
computation  of  inertia  properties  for  RBI  and  RB2  followed  the  same  procedure  as 
the  wings.  Table  5.2  shows  the  hnal  dimensionless  inertia  properties  for  the  model 
components. 

In  order  to  complete  the  set  of  parameters  required  for  the  FSI  simulation  of 
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longitudinal  flight  the  dimensionless  held  acceleration  (gravity)  was  set  to  g  =  0.01, 
a  value  that  gives  a  weight  force  lower  than  the  required  to  balance  the  mean  lift 
force  of  the  prescribed  simulation.  Thus,  it  is  expected  that  the  model  would  have 
an  upward  net  motion.  Also  the  initial  conditions  on  the  state  variables  were  set 
to  Xo  =  0,  Zo  =  —0.56,  9o  =  90°,  620  =  0,  and  zero  velocities.  A  unit  dimensionless 
torsion  stiffness  Kt  =  1  was  employed  at  the  abdomen  hinge.  This  was  seen  to 
maintain  the  relative  angle  between  RBI  and  RB2  within  a  5"  amplitude.  An  initial 
attempt  using  the  kinematics  of  the  previous  prescribed  motion  simulation,  where 
the  mean  stroke  angle  is  5°  behind  the  b2  axis,  gave  a  large  pitch  down  acceleration, 
due  to  the  moment  imbalance  with  respect  to  y.  In  the  simulation  we  present  below, 
we  changed  the  mean  stroke  angle  to  'ipm  =  5.57r/180,  advanced  respect  to  the  b2 
axis.  In  this  case  a  slight  pitch  up  motion  was  found  along  the  simulation.  It  is 
important  to  note  that,  in  several  flapping  wing  systems,  both  hovering  and  forward 
flight  have  been  found  to  be  dynamically  unstable  through  linear  stability  analysis 
(i.e.  [94],  [97]).  Given  the  system  of  aerodynamic  forces  produced  through  flapping, 
and  regardless  how  balanced  the  mean  force  and  moment  equations  are  initially,  it 
is  expected  that  the  solution  will  eventually  diverge  as  the  integration  progresses. 

In  hgure  5.17  the  values  of  the  state  variables  x(t),  y(t),  6{t)  and  6*2(t)  as  a 
function  of  integration  time  are  shown.  The  vertical  position  z{t)  is  seen  to  increase 
throughout  the  calculation  consistent  with  the  fact  that  the  mean  resulting  vertical 
force  is  directed  upwards.  The  horizontal  coordinate  of  the  center  of  mass  of  RBI 
takes  oscillating  positive  values,  which  diminish  as  the  calculation  progresses.  This 
is  due  to  the  fact  that  the  orientation  angle  d{t)  of  RBI  increases  with  time  and  so 
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does  the  stroke  plane  angle  6{t)  +  /3-,  introducing  a  component  of  the  force  normal 
to  the  stroke  plane  angle  (lift  in  the  prescribed  simulation)  along  the  negative  x 
direction.  The  oscillatory  component  of  x{t)  is  about  0.16,  which  is  consistent  to 
the  data  reported  in  [112]  for  a  Manduca  Sexta  hawkmoth  model  with  similar  wing 
to  body  mass  ratio.  The  orientation  angle  6{t)  increases  steadily,  due  to  moment 
imbalance.  A  lower  value  of  'ipm  is  required  to  reduce  this  effect.  Also,  9(t)  oscillates 
with  a  peak  to  peak  amplitude  of  4°  similarly  to  what  is  reported  in  reference  [112]. 


0  0.5  1  1.5  2  2.5  3  3.5 

t/T 


Figure  5.17:  Variation  of  positions  x(t),  z(t)  and  angular  coordinates  9(t)  (blue), 
02 (t)  (red)  with  time  for  FSI  simulation  of  free  longitudinal  flight. 

In  hgures  5.18  the  positions  of  the  center  of  mass  of  RBI  are  shown,  together 
with  the  location  of  the  thorax,  abdomen  and  the  orientation  of  the  stroke  plane  at 
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different  time  points.  The  increase  of  the  angle  6{t)  as  the  simulation  progresses  is 
given  by  a  pitch  up  momentum  of  the  force  normal  to  the  stroke  plane.  Aerodynamic 
forces  are  follower  forces  in  the  sense  that  they  are  always  aligned  (i.e.  are  normal 
or  tangent)  to  the  aerodynamic  surfaces,  and  their  dynamic  effect  is  in  general 
destabilizing.  This  fact  points  to  the  need  of  developing  control  strategies  using 
some  of  the  wing  kinematics  parameters  described  (i.e.  'ipm,  «)  in  order  to  minimize 


pitch  oscillations  or  divergence. 
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Figure  5.18:  Variation  of  position  of  the  center  of  mass  of  RBI  with  time.  The  model 
section  on  the  plane  fii  —  hs  is  shown  along  with  the  orientation  of  the  stroke  plane 
(in  red)  for:  (a)  t/T  =  0.73,  (b)  t/T  =  1.38,  (c)  t/T  =  2.13,  and  (d)  t/T  =  2.80. 
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Figure  5.19:  Variation  of  force  coefficients  C^it),  and  C'Ar(t)  (force  normal 

to  stroke  plane)  with  time  for  FSI  calculation  of  flapping  wing  model  in  free  longi¬ 
tudinal  flight:  (blue)  RBI,  (dashed  magenta)  RB2,  (red)  RW,  and  (•)  LW.  The 
green  curves  correspond  to  C^it),  and  CL{t)  from  the  prescribed  kinematics 

simulation. 

In  figures  5.19  the  time  variation  of  force  coefficients  in  the  horizontal  x,  and 
vertical  z  directions,  and  also  in  a  direction  normal  to  the  stroke  plane  (lift  direction 
for  the  prescribed  kinematics  simulation)  are  shown.  We  see  that  the  fluid  forces 
computed  for  RW  and  LW  are  again  equal.  Also,  as  in  the  prescribed  kinematics 
simulation  the  force  in  the  y  direction  was  found  to  be  negligible  with  respect  to 
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the  other  components.  The  maximum  force  on  the  x  direction  provided  by  RBI 
and  RB2  is  8%  the  force  on  the  wings.  It  is  seen  that,  due  to  the  stroke  plane 
rotation,  Cx  and  Cz  vary  signihcantly  as  with  time.  The  coefficients  resulting  from 
the  prescribed  kinematics  calculation  are  plotted  on  the  same  hgures  for  comparison. 
All  force  coefficients  take  lower  values  than  the  hxed  body  calculation.  This  is  due 
to  the  fact  that  body  motion  lowers  the  ability  of  the  wings  to  transfer  momentum 
to  the  fluid  (see  also  [112]). 

The  instantaneous  flow-helds  at  different  simulation  times  are  seen  in  hg¬ 
ures  5.20.  Here,  an  isocontour  of  Q  colored  by  the  vorticity  Uy  shows  the  evolution 
of  the  fundamental  how  structures,  namely  leading  edge  and  tip  vortices.  In  hgures 
5.20a-b  we  see  the  leading  edge  vortices  attached  to  each  wing,  which  detach  in 
the  vicinity  of  the  wing  tips.  Here  vorticity  is  reoriented  forming  the  wing  tip  vor¬ 
tices.  Vorticity  is  also  shed  from  the  regions  of  the  wing  proximal  to  the  bodies,  and 
therefore,  two  vortical  structures  are  being  generated  on  each  wing  (hgures  5.20c-d). 
This  secondary  vortices  are  dependent  on  the  planform  geometry  of  the  wings  used. 
The  vertical  displacement  of  the  model  due  to  the  imbalance  in  the  lift  force  and 
its  own  weight  can  be  clearly  seen  in  the  sequence  of  instantaneous  snapshots  in 
hgures  5.20. 
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Figure  5.20:  Q  isocontour  colored  by  vorticity  on  the  y  direction.  40  contours  of  ujy 
from  -20  to  20  are  used,  (a)-(h):  t/T  =  0.75,  1,  1.25,  1.5,  1.75,  2,  2.25  and  2.5. 
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Chapter  6 


Summary,  Contributions  and  Directions  for  Future  Work 

In  this  work,  an  adaptive  mesh  refinement,  large  eddy  simulation  strategy  has 
been  developed  for  fluid-structure  interaction  problems  in  transitional  and  turbu¬ 
lent  flow  regimes.  An  immersed  boundary  reconstruction  procedure  to  represent  the 
moving  and/or  deforming  bodies  immersed  in  the  flow  has  been  proposed.  The  over¬ 
all  method  is  a  generalization  of  the  formulation  initially  proposed  by  Uhlmann  [101], 
where  the  main  difference  with  existing  direct-forcing  schemes  (i.e.  [37],  [50],  [9]) 
is  that  the  evaluation  of  the  forcing  function  is  done  on  the  Lagrangian  markers 
instead  of  the  Eulerian  points.  The  main  advantages  of  this  strategy  compared  to 
existing  direct-forcing  schemes  can  be  summarized  as  follows: 

i)  It  is  more  versatile  since  it  decouples  the  local  discretization  from  the  compu¬ 
tation  of  the  forcing  function  and,  therefore,  can  be  implemented  into  struc¬ 
tured  or  unstructured  codes  in  a  straightforward  manner.  Most  of  the  available 
direct-forcing  schemes  have  been  developed  in  the  framework  of  finite-difference 
or  finite-volume  formulations  on  Cartesian  grids.  The  proposed  scheme  can  be 
used  with  other  spatial  discretization  approaches  (i.e.,  finite  elements). 

ii)  It  is  very  robust  in  dealing  with  contact  and  collisions  of  multiple  bodies.  The 
forcing  function  is  built  and  appropriately  scaled  based  on  the  contributions  of 
all  bodies  in  the  vicinity  of  an  Eulerian  point  without  any  special  treatment.  In 
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most  direct-forcing  methods,  the  presence  of  two  or  more  Lagrangian  markers 
from  different  bodies  in  the  proximity  of  the  same  Eulerian  grid  point  is  usually 
the  source  of  ambiguity. 

The  immersed  boundaqry  method  was  also  found  to  maintain  the  second-order  spa¬ 
tial  accuracy  of  the  underlying  hnite-difference  solver.  Most  importantly,  it  has  been 
demonstrated  that  when  combined  with  the  scheme  that  has  been  proposed  to  com¬ 
pute  the  surface  forces,  it  has  a  sharp-like  behavior  similar  to  sharp  Eulerian  direct- 
forcing  schemes  and  boundary  conforming  methods.  The  overall  computational  cost 
of  the  proposed  formulation  is  comparable  to  other  direct  forcing  approaches  avail¬ 
able  in  the  literature  (i.e.  [9],  [50]).  In  the  case  of  stationary  bodies,  the  interface 
tracking  step,  as  well  as  the  computations  of  the  shape  functions  can  be  carried  out 
in  a  pre-processing  module  and  stored  in  memory.  Then,  the  computational  effort 
associated  with  the  forcing  step  is  reduced  to  a  small  fraction  of  the  overall  cost. 
In  the  case  of  moving  bodies  the  operations  described  in  Section  2.3.1  need  to  be 
performed  at  each  timestep.  The  overall  cost  depends  on  the  total  number  of  La¬ 
grangian  markers  and  the  size  of  the  Eulerian  grid.  The  tracking  step  in  the  current 
formulation,  is  probably  less  expensive  than  other  direct-forcing  schemes  since  only 
the  closest  point  to  each  marker  needs  to  be  identihed.  The  forcing  step  on  the 
other  hand,  is  more  expensive  than  typical  direct-forcing  schemes  since  the  shape 
function  computation  involves  the  solution  of  a  4  x  4  linear  system  (equation  2.11). 
However,  like  all  other  direct-forcing  schemes,  the  cost  per  Lagrangian  marker  is 
constant.  Therefore,  the  overall  cost  is  proportional  to  the  number  of  Lagrangian 
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markers. 


To  extend  the  range  of  applicability  of  immersed  boundary  methods  such  as 
the  one  proposed  above,  to  high  Reynolds  number  complex  flows,  the  author  has 
also  developed  a  structured  AMR  method.  In  this  approach,  a  single-block  solver  is 
employed  on  a  hierarchy  of  sub-grids  with  varying  spatial  resolutions.  Each  of  these 
sub-grid  blocks  has  a  structured  Cartesian  topology,  and  each  block  is  part  of  a  tree 
data  structure  that  covers  the  entire  computational  domain.  Time  advancement  is 
done  by  using  a  fractional  step  method.  All  spatial  derivatives  are  approximated 
with  second  order  hnite-differences  on  a  staggered  grid.  The  Paramesh  toolkit  [62] 
is  utilized  to  keep  track  of  the  grid  hierarchy,  and  perform  the  required  restric¬ 
tion/prolongation  and  guard-cell  hlling  operations.  The  author  has  demonstrated 
that  the  accuracy  of  dynamic  AMR  is  greatly  affected  by  the  conservation  prop¬ 
erties  of  the  prolongation  and  restriction  operators.  The  author  has  developed  a 
divergence-preserving  prolongation  operator  tailored  to  the  specihc  AMR  topology, 
where  the  grid  size  between  consecutive  rehnement  levels  can  only  differ  by  a  fac¬ 
tor  of  two.  Overall,  the  second-order  spatial  and  temporal  accuracy  of  the  basic 
solver  are  maintained.  The  computational  efficiency  of  the  proposed  formulation  is 
a  balancing  act  between  the  lower  CPU  and  memory  cost  incurred  bn  using  AMR, 
and  the  additional  work  derived  from  the  augmented  computational  complexity.  In 
general,  the  proposed  solver  has  a  signihcant  advantage  over  single  block  solvers,  in 
high  Reynolds  number  FSI  problems,  where  the  hne  grid  patches  needed  to  resolve 
the  boundary  layers  on  the  body  surfaces  have  to  be  constantly  rearranged  follow¬ 
ing  the  motion.  For  the  falling  plates  simulation  of  Section  3.2.3,  for  example,  the 
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number  of  grid  points  utilized  in  the  AMR  computation  is  1/10  of  the  ones  used  in 
the  uniform  grid  solution.  In  terms  of  wall  time,  the  AMR  calculation  took  about 
3  seconds  per  timestep  on  a  single  processor,  while  the  equivalent  uniform  grid  run 
required  roughly  twice  that  amount.  It  is  also  noted  that  the  uniform  grid  solver 
utilizes  a  direct  solver  for  the  Poisson  equation  (i.e.  [9]),  which  is  much  more  effi¬ 
cient  compared  to  the  multigrid  solver.  If  a  multigrid  solver  is  to  be  adopted  for  the 
case  of  the  uniform  grid  too,  then,  the  computational  savings  by  using  AMR  would 
have  been  an  order  of  magnitude  higher. 

In  general,  it  is  found  for  projection  schemes  for  the  Navier-Stokes  equations 
that  the  Poisson  equation  solution  takes  a  substantial  part  of  the  computing  cost. 
The  problem  is  particularly  acute  for  parallel  multigrid  solution  schemes  on  octree 
meshes,  where  sequential  relaxation  across  rehnement  levels  requires  extensive  data 
communication  among  processors.  Also,  work  load  balancing  among  processors 
becomes  a  very  difficult  task.  The  use  of  a  hybrid  multigrid  and  a  highly  parallel 
direct  solver  for  the  coarse  grid  solution,  together  with  sizable  meshes  on  the  coarsest 
level,  is  found  to  improve  the  computation  times  on  some  problems,  for  example,  the 
computations  of  flow  around  a  sphere.  It  is  believed  that  highly  scalable  Poisson 
solvers  for  AMR  grid  structures  like  the  one  employed  in  this  work,  need  to  be 
developed  to  perform  the  next  generation  of  numerical  simulations. 

The  author  has  also  shown  that  the  proposed  AMR  scheme  is  well  suited  for 
carrying  out  LES  studies.  In  Chapter  4,  we  performed  numerical  simulations  of 
spatially  developing  homogeneous  isotropic  turbulence,  convected  through  an  in¬ 
terface  where  the  grid  is  suddenly  coarsened  or  rehned  by  a  factor  of  two  in  each 
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direction.  The  author  has  compared  calculations  where  the  hlter  width  is  propor¬ 
tional  to  the  grid  size  and  is  discontinuous  at  the  coarse-hne  grid  interface,  with 
cases  in  which  the  hlter-width  varies  gradually  between  the  values  corresponding 
to  the  coarse  and  hne  grids.  Both  the  Smagorinsky  and  the  Lagrangian-Dynamic 
Eddy  viscosity  (LDEV)  SGS  models  have  been  used.  In  addition,  the  explicit  hlter- 
ing  of  the  advective  term  in  conjunction  with  the  application  of  the  LDEV  model 
has  been  evaluated.  A  sudden  refinement  of  the  grid  does  not  result  in  a  significant 
flow  perturbation;  small  scales  are  gradually  generated  downstream  of  the  interface, 
resulting  in  a  smooth  flow  across  the  interface.  When  the  grid  is  suddenly  coars¬ 
ened,  on  the  other  hand,  a  considerable  energy  pileup  at  small  scales  is  observed 
near  the  interface.  For  coarse-hne  interfaces  a  discontinuous  hlter  width  gives  more 
accurate  results  than  a  smoothly  decreasing  one:  decreasing  the  eddy  viscosity  al¬ 
lows  small  eddies  to  be  generated  more  rapidly,  while  the  smoothly  varying  hlter 
gives  increased  eddy  viscosity  on  the  hue-grid  side,  which  delays  the  generation  of 
small  scales.  When  the  grid  is  suddenly  coarsened,  on  the  other  hand,  an  increase 
of  the  eddy  viscosity  upstream  of  the  interface  through  a  smooth  increase  of  the 
hlter-width  is  found  to  be  benehcial.  For  hne-to-coarse  interfaces,  the  LDEV  model 
with  a  smoothly  varying  interface  is  found  to  give  more  accurate  results  than  the 
Smagorinsky  model,  and  the  how  transitions  to  the  single-grid  results  within  one 
integral  scale  Ln.  The  interface  ehect  is  inherently  related  to  the  amount  of  eddy 
viscosity  modeling  employed  in  the  simulations.  As  the  Reynolds  number  Re\j  is 
increased,  for  the  hne-coarse  transition,  the  intensity  and  extension  of  the  energy 
concentration  zone  on  the  hne  side  is  magnihed.  Correspondingly,  on  the  coarse-hne 
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situation,  the  distance  required  for  the  restoration  of  equilibrium  is  enhanced  as  the 
resolution  becomes  coarser. 

Explicit  hltering  has  the  the  ability  to  decrease  the  high  wavenumber  energy 
content  of  the  flow  prior  to  the  interface,  and  this  improves  the  transition  across  the 
grid  discontinnities.  This  strategy  works  very  well  with  the  present  AMR  scheme, 
and  was  fonnd  to  give  excellent  resnlts  for  the  case  of  the  flow  aronnd  a  sphere  at 
Re  =  10^.  The  fact  that  the  LDEV  model,  which  has  some  memory  effects,  gives 
better  results  than  the  Smagorinsky  one  indicates  that,  perhaps  models  that  inclnde 
a  transport  equation  (which  also  would  include  memory  effects)  may  be  benehcial. 
Further  investigation  of  other  LES  modeling  strategies  for  AMR  is  left  as  a  poten¬ 
tial  fntnre  work  direction.  Compared  to  available  strnctnred  [64]  or  nnstrnctnred 
AMR  [45]  formulations,  the  constraint  that  neighboring  blocks  can  only  differ  by 
one  level  of  refinement  may  result  in  larger  overall  grids  to  achieve  the  same  local 
resolntion,  especially  in  internal  flows.  This  is  not  a  major  concern  in  LES  studies 
of  turbulent  and  transitional  flows,  where  grid  discontinnities  generated  by  neigh¬ 
boring  blocks  that  differ  by  more  than  a  factor  of  two  can  contaminate  the  high 
freqnency  content  and  nnphysically  enhance  tnrbnlent  flnctnations  on  the  hne  grid 
side;  therefore,  they  are  not  desirable. 

In  Chapter  5,  the  nnmerical  tools  were  applied  to  study  and  simnlate  prob¬ 
lems  relevant  to  flapping  wing  systems  at  low  Reynolds  numbers  regimes.  First,  the 
influence  of  flexibility  on  the  aerodynamic  performance  of  a  two-dimensional  hov¬ 
ering  wing  section  was  nnmerically  stndied.  Here,  the  wing  model  consists  of  two 
rigid  links  that  are  joined  at  the  center  with  a  linear  torsion  spring.  By  prescribing 
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the  kinematics  of  the  top  link,  the  structural  system  is  effectively  reduced  to  a  sin¬ 
gle  degree-of-freedom  non-linear  oscillator.  The  results  obtained  demonstrate  that 
flexibility  can  be  benehcial  in  terms  of  enhancing  the  aerodynamic  performance. 
Furthermore,  it  has  been  found  that  in  the  frequency  range  below  the  hrst  natural 
frequency,  the  best  performance  is  achieved  when  the  wing  is  driven  at  a  frequency 
close  to  one  of  the  non-linear  resonances  (a  superharmonic  resonance  of  order  three) 
of  the  system.  This  behavior  is  common  to  all  of  the  Reynolds  numbers  studied.  In 
terms  of  the  flow  physics,  the  wake  capture  mechanism  is  enhanced  partially  due  to 
a  stronger  flow  around  the  wing  at  stroke  reversal.  However,  it  needs  to  be  noted 
that  the  cases  where  the  wing  is  driven  at  or  close  to  the  hrst  natural  frequency 
of  the  system  were  not  considered  in  this  study,  and  it  is  possible  that  a  better 
aerodynamic  performance  may  be  achieved  at  the  fundamental  resonance  and  this 
remains  to  be  explored.  The  study  also  leads  to  the  following  open  questions: 

i)  Why  is  there  a  performance  enhancement  when  the  system  is  excited  at  a 
happing  wings  non-linear  resonance  and  would  one  achieve  a  better  performance 
with  a  non-linear  resonance  compared  with  a  linear  resonance? 

ii)  Which  kinematics  is  preferable  from  an  aerodynamic  efficiency  standpoint? 

The  interplay  between  wing  hexibility  and  kinematics  together  with  qualitative 
changes  in  the  system  dynamics  as  a  function  of  the  Re  number  requires  further 
investigation,  and  this  is  left  as  a  direction  for  future  research. 

A  four  rigid  body,  happing  wing  model  has  also  been  developed.  The  model  is 
composed  of  the  thorax-head,  abdomen,  and  a  pair  of  wings.  The  kinematics  for  the 
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general  motion  in  three  dimensions  has  been  studied,  and  the  equations  of  motion 
have  been  derived  for  the  simple  case  of  free  longitudinal  flight.  A  simple  strategy  to 
estimate  the  inertia  properties  of  each  of  the  components  has  been  developed.  Both 
all  prescribed  kinematics  and  FSI  simulations  of  longitudinal  flight  were  performed, 
showing  that  the  proposed  numerical  methodology  is  well  suited  to  approach  these 
complex  problems.  It  is  believed  that  this  work  will  serve  as  a  basis  for  more  in-depth 
studies  into  dynamics,  stability  and  control  in  free  flight  regimes. 
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Appendix  A 


Divergence  preserving  prolongation  in  three  dimensions 

The  development  of  divergence  preserving  operators  in  three-dimensions  is 
similar  to  the  one  for  two-dimensions  outlined  in  section  3.1.2.  In  particular,  the 
discrete  divergence  of  the  coarse-grid  cell  shown  in  Fig.  A.l  is: 


To  interpolate  the  velocities  on  the  cell  faces,  the  one-dimensional  interpolations 
in  two-dimensions,  are  now  replaced  with  two-dimensional  quadratic  interpolations 
assuming  a  polynomial  variation  of  the  form; 


7])  =  ao  +  ai^  -h  a2?7  -h  as^r]  +  04^^  -h  a^rf  +  (A. 2) 


An  example  stencil  is  shown  in  hgure  A. la,  and  a  system  analogous  to  (3.4)  can 
be  constructed  by  applying  Eq.  (A. 2)  in  the  o:,  jS  =  —1,0,1  positions 

at  level  /,  except  for  {^i,'r]j)-  Then  the  known  value  of  the  coarse  grid  variable  at 
location  {^i,rij)  is  used  in 


7/^  =  —  -I- 7;^~*~^  -I- 7;^~*~^ 


+  1,V-|-1 


-L  7/^+1 


(A.3) 


and  the  four  hne  grid  values  are  computed  by  interpolations  from  Eq.  (A. 2).  The  ma¬ 
trix  associated  to  the  resulting  system  can  be  inverted  using  symbolic  manipulation 
software,  such  that  the  computation  of  the  interpolation  coefficients  a*,  i  =  0, 1, . . . ,  8 
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(b) 


Figure  A.l;  (a)  Interpolation  stencil  for  2D  face  interpolations  used  in  three  di¬ 
mensional  divergence  preserving  prolongation,  (b)  Internal  variables  to  a  given 
coarse-grid  cell  indexed  by  i,  j,  k. 


will  require  a  9  x  9  matrix-vector  multiplication.  The  computation  of  the  four  fine 
grid  variables  requires  four  additional  9x1  vector-vector  multiplications. 

Once  the  velocity  field  has  been  found  on  the  cell  faces,  the  velocities  in¬ 


ternal  to  the  coarse  cell  need  to  be  obtained.  There  are  twelve  internal  hne-grid 
vaiiauies,  u-j/_ijv_|_i^fc/,  «j/_i 

that  need  to  be  determined 

(see  Fig.  A.lb).  Four  of  the  twelve  4-\,y,fco  4tAfc'+i)’ 

can  be  found  using  one  dimensional  quadratic  interpolation  from  the  corresponding 


coarse-face  variables,  and  for  the  remaining  ones,  eight  discrete  divergence  equations 
are  used  as  in  the  two-dimensional  example.  The  resulting  system  can  be  written 
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as: 


hand  side  are  given  by 


,  _  ^  _  ^i'-l,j',k'  ^i'-2,j',k'  _  ^i'-l,j',k'  ^i'-l,j'-l,k'  ,k'-l 

^^+1  _  ^^+1 

7  r-i  ^i',j',k'  “i'-lj'.fc'  ,  ^i',j'-l,k'  ,  ^i',j',k'-l 

^2  =  - 


A^^+1 


— 

Ax^+i  A|/^+^  Az^+i 

^,^  +  1  _  7/  +  I 


„.(+!  „y+i  _  „y+i 

,  r-i  ,  “7'-2,?''+1,A:'  ^i' -l,j' +l,k'  ^i'-l,j', 

'>3  =  - - 


A|/^+i 


'^i'-ij'+i,fc'-i 

A2:'+^ 


^  ^  ^  Az'+i 


^  Aa;'+^  A|/'+^ 

^,^+1  „1^+l  _  „1^  +  l  r,,J-  +  ^ 


^y^+1  y/+l  _  y/  +  l 

Oe  =  -Du - -^2^— - - ir-fTT — ^ - ^ 

Ax'+i  A|/'+i  Az 


A2:^+^ 


+  - £ 


L 

1  o'-Ll  Jh/_U1 


(A.5) 


&8 


n  “i'j'+l.fc'+l  “i'-lJ'+l,A:'+l  A',/+l,fc'+l  ‘^i',j',fc'+l 

^  A^DI 

7,/+! 

^i'j'+l,fc'+l 

A2:*+^ 
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As  this  system  is  of  rank  7,  one  more  equation  similar  to  (3.6)  is  added.  As  in 
the  two-dimensional  case,  the  left-pseudo-inverse  of  the  resulting  9x8  coefficient 
matrix  Ae  can  be  computed  symbolically.  The  implementation  requires  additional 
8x9  matrix-vector  operations  to  obtain  the  unknown  variables.  It  is  possible  to 
reduce  the  size  of  system  (A. 4),  by  replacing  the  values  for  and 

given  by  the  hrst  and  last  equations,  on  the  remainder  ones.  Then,  the  condensed 
system  is  of  size  6x6  and  still  requires  an  extra  equation  to  be  full  rank.  In 
this  case  the  computer  implementation  would  require  a  smaller  6x7  matrix-vector 
multiplication  to  obtain  the  rest  of  the  hne-grid  variables. 
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